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Abstract 


This  research  includes  an  investigation  of  the  mechanisms  of  diffraction  and 
reinitiation  that  enable  a  detonation  diffuser.  It  describes  a  set  of  geometric  parameters 
necessary  to  design  a  diffuser  for  a  given  detonable  mixture  and  initial  channel  height. 
Predetonators  with  channel  height  less  than  the  critical  height  are  ineffective  because 
detonations  in  small  channels  decouple  into  separate  shock  and  combustion  fronts  when 
the  channel  height  increases.  A  detonation  diffuser  allows  the  channel  height  to  increase 
by  utilizing  the  decoupled  shock  wave  to  reinitiate  detonation.  In  the  diffuser,  a 
detonation  initially  decouples  into  separate  shock  and  combustion  fronts,  and  then  the 
decoupled  shock  front  reflects  from  an  oblique  surface  initiating  a  secondary  detonation 
that  survives  the  expansion.  This  research  investigated  the  three  regions  of  a  detonation 
diffuser:  the  initial  diffraction,  the  reflecting  surface,  and  the  second  diffraction  comer. 
Schlieren  video  of  two-dimensional  diffracting  detonations  recorded  the  position  of  the 
detonation,  decoupled  shock  front  and  flame  front.  Observations  of  the  decoupled  shocks 
reflecting  from  surfaces  showed  that  a  45°  reflecting  surface  must  be  placed  less  than  80 
mm  downstream  of  the  initial  diffraction  comer  to  initiate  a  secondary  detonation  in  more 
than  91%  of  repeated  trials.  Observations  of  the  interaction  of  diffracting  detonations 
with  multiple  obstacles  revealed  that  the  best  performance  (smallest  separation,  and 
highest  Mach  number)  occurred  when  the  decoupled  shock  reflected  from  four  separate 
obstacles  at  approximately  the  same  time. 
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I.  Introduction 

Motivation 

Failure  of  a  detonation  front  in  the  transition  from  a  subcritical  channel  to  a  supercritical 
channel  is  an  unaddressed  area  of  concern  in  the  design  of  pulsed  detonation  engines.  When  a 
predetonator  is  used  to  initiate  a  detonation,  the  predetonator  channel  should  be  as  small  as 
possible  to  minimize  the  requirements  for  reactants,  but  the  thrust  tube  of  the  PDE  should  be  as 
large  as  possible  to  maximize  the  thrust  per  cycle.  The  detonation  in  the  predetonator  channel 
fails  when  the  shock  and  combustion  fronts  decouple  due  to  the  area  increase.  The  once 
decoupled,  the  combustion  is  less  efficient  than  a  detonation  and  raises  the  pressure  less  than  the 
desired  detonation.  The  benefits  of  detonation  will  be  maintained  if  the  detonation  is  reinitiated 
after  decoupling.  A  detonation  diffuser  is  a  device  designed  to  reinitiated  detonation  during  the 
transition  from  a  subcritical  channel,  such  as  the  predetonator,  to  a  supercritical  channel,  such  as 
the  thrust  tube.  The  detonation  diffuser  will  utilize  the  decoupled  shock  front  to  reinitiate 
detonation  as  the  height  of  the  channel  increases. 

At  a  sudden  area  expansion,  diffracting  detonations  decouple  or  not  depending  on  the 
initial  channel  height  (Zeldovich,  1956).  The  critical  channel  height,  which  depends  on  the  cell 
size  of  the  reactant  mixture,  determines  whether  decoupling  occurs  (Mitrofanov,  1965).  In 
supercritical  channels,  the  initial  channel  height  is  greater  than  the  critical  height,  and  the 
detonation  diffracts  without  decoupling  (Fig.  la).  In  subcritical  channels,  the  detonation 
decouples  into  separate  shock  and  deflagration  fronts  (Fig.  Ic).  At  the  critical  height,  decoupling 
occurs  initially,  but  naturally  occurring,  localized  explosions  reinitiate  detonation  in  the  space 
between  the  shock  and  deflagration  fronts  and  restore  the  detonation  mode  of  combustion 
(Soloukhin  and  Ragland,  1969)  (Fig.  lb). 
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a)  Supercritical:  b)  Critical:  c)  Subcritical: 

No  decoupling/  Decoupling  followed  Complete  decoupling/ 

successful  by  reinitiation  failed  detonation 


Figure  1.  Diffraction  regimes  in  sudden  area  expansion 
Methods  will  be  discussed  later  that  involve  using  shock  reflections  to  restore  a 

diffracting  detonation.  Shock  reflection  causes  a  local  explosion  by  compressing  the  reactants  so 
rapidly  that  chemical  reaction  begins  and  remains  coupled  to  the  reflected  shock  wave  (Brown 
and  Thomas,  2000).  Brown  and  Thomas  (2000)  suggested  that  the  presence  of  a  boundary  layer 
is  necessary  for  shock  initiation,  but  Thomas  et  al.  (2002)  demonstrated  that  shock  compression 
alone  is  sufficient  to  cause  localized  explosions  and  initiate  detonation  by  reflecting  a  non¬ 
reacting  shock  with  the  end  of  a  cylinder  (Fig.  2).  The  cylinder  experiment  eliminated  the 
boundary  layer  interactions  that  were  present  in  Brown  and  Thomas  (2000).  The  minimum 
shock  strength  for  localized  explosion  depends  on  the  reflecting  surface  area,  speed  of  sound  in 
the  undisturbed  reactants,  and  ignition  delay  (Thomas  et  al.,  2002).  Thomas  et  al.  defined  a 
criterion  for  detonation  initiation  based  on  these  properties  that  will  be  discussed  in  Chapter  II. 
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Figure  2.  Detonation  initiated  by  non-reacting  shock  (Thomas  et  al.,  2002) 

The  minimum  incident  Mach  number  for  initiation,  reported  by  Thomas  et  al.,  is  2.7.  The 

Chapman-Jouguet  Mach  number  is  5.3,  and  there  is  sufficient  potential  in  the  shockwave  from  a 
recently  decoupled  detonation  to  initiate  a  new  detonation.  Reinitiation  can  be  achieved  by 
reflecting  the  decoupled  shock  a  done  by  Thomas  et  al.  did.  From  this  reasoning,  it  seemed 
possible  to  construct  a  detonation  diffuser  utilizing  reinitiation  of  the  decoupled  detonations. 
Because  the  diffuser  requires  only  a  sufficiently  strong  shock  to  initiate  detonation  it  functions 
even  when  the  initial  channel  height  is  subcritical.  Unlike  the  cylinder  in  Figure  2,  a  detonation 
diffuser  must  reinitiate  detonation  in  the  limited  time  between  the  passing  of  the  decoupled  shock 
and  combustion  fronts.  Normal  reflection  of  the  shock  would  result  in  a  detonation  that  runs  out 
of  reactants  when  it  encounters  the  combustion  front.  Rotating  the  reflecting  surface  such  that 
the  shock  reflection  is  oblique  preserves  most  of  the  compression  gained  by  the  reflection  while 
giving  the  newly  formed  detonation  front  a  route  to  escape  the  oncoming  combustion  front.  This 
research  investigates  the  reflecting  angle  and  position  relative  to  a  decoupling  detonation  of 
reflecting  surfaces.  The  goal  is  to  induce  a  planned  localized  explosion  and  reinitiate  a 
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detonation  that  decoupled  due  to  sub-critical  diffraction.  As  will  be  seen,  the  initial  diffraction 
decouples  the  shock  and  combustion  in  the  primary  detonation.  The  decoupled  shock  reflects 
from  a  reflecting  surface  causing  a  localized  explosion  that  evolves  into  a  secondary  detonation 
(see  Fig.  3). 

The  secondary  detonation  in  turn  diffracts  at  a  second  step  and  reinitiates  at  a  second 
reflecting  surface.  The  cycle  of  diffraction  and  reinitiation  continues  until  the  channel  height 
exceeds  the  critical  height.  Stevens  et  al.  (2011)  published  the  first  example  of  a  single  step 
detonation  diffuser  (Fig.  3).  The  diffuser  employed  a  converging  ramp  as  the  reflecting  surface. 
The  ramp  angle  (P)  was  14°,  with  a  vertical  offset  of  13  mm  and  a  rise  of  13  mm  resulting  in 
zero  net  expansion.  Local  explosions  occurred  near  the  middle  of  the  converging  ramp  where 
the  expansion  ratio  (Eq.  1)  was  1.17. 


Expansion  ratio  = 


Channel  height  at  local  explosion 
Initial  channel  height 


(1) 


Figure  3.  Converging  ramp  configuration  used  by  Stevens  et  al.  (2011) 

A  detonation  diffuser  such  as  one  shown  in  Fig.  3  is  applicable  to  any  situation  requiring 

detonation  in  a  supercritical  channel  that  is  supplied  by  a  subcritical  channel  such  as  the 
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transition  between  predetonator  and  thrust  tube  in  PDEs.  Ideally,  the  predetonator  channel  is 
subcritical  to  minimize  the  volume  of  sensitive  mixture,  and  the  thrust  tube  is  as  large  as  possible 
to  maximize  the  thrust  per  pulse.  The  predetonator  is  a  promising  initiation  means  due  to  small 
volume  and  extremely  short  detonation  initiation  times  and  distance,  but  predetonators  see 
limited  use  due  to  sporadic  transmission  of  the  pre-detonator  detonation  to  the  thrust  tubes 
(Hoke,  2006).  A  detonation  diffuser  will  remedy  the  sporadic  transmission. 

Research  Objectives 

This  research  investigates  the  feasibility  of  a  detonation  diffuser.  The  results  demonstrate 
and  parameterize  direct  initiation  of  a  secondary  detonation  by  the  reflection  of  a  shockwave 
formed  when  a  detonation  decouples  at  the  exit  of  a  subcritical  channel.  The  diffuser  design  is 
built  upon  the  design  studied  by  Stevens  et  al.  (2012)  and  shown  in  Figure  3.  This  research 
studies  the  initial  decoupling  of  the  detonation  exiting  the  subcritical  channel  to  determine  the 
locations  where  the  shock  propagation  Mach  number  is  sufficient  for  reinitiation  and  the 
locations  where  the  decoupled  flame  front  prevents  the  secondary  detonation  from  propagating  to 
the  exit  of  the  diffuser.  This  research  also  studies  the  initiation  of  secondary  detonations  to 
determine  the  reflecting  surface  angle  and  position  that  result  in  initiation  of  a  secondary 
detonation.  Finally,  this  research  studies  a  series  of  diffuser  configurations  to  determine  what 
effect  the  number  of  reflecting  surfaces  and  their  arrangement  has  on  the  formation  and  survival 
of  secondary  detonations. 

This  research  examines  initial  diffraction,  reinitiation  (initiation  of  secondary 
detonations),  and  secondary  detonation  decoupling  in  turn  to  identify  and  bound  the  important 
parameters.  Figure  4  shows  the  regions  of  interest  for  initial  diffraction,  reinitiation,  and 
secondary  detonation  propagation.  In  the  first  phase,  the  initial  diffraction  at  the  first  diffraction 
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comer  reveals  the  decoupled  shoek  strength  and  distanee  between  the  shock  and  combustion 
fronts.  Knowledge  of  the  shock  Mach  number  (Mshock  =  Vshock/a)  is  neeessary  to  position  and 
orient  the  reflecting  surface  such  the  initiation  eriterion  is  satisfied. 


- Phase  1 .  Initial  Diffraetion 

- Phase  2.  Reinitation 

- Phase  3:  Secondary  Detonation 

_ Propagation _ 


First  Reflecting  Second 


Figure  4.  Sequential  order  and  general  position  of  phenomena  in  a  detonation  diffuser 

Knowledge  of  the  separation  distanee  between  shoek  and  eombustion  or  “shock-flame 
separation”  is  neeessary  to  avoid  trapping  the  seeondary  detonation  between  the  first  diffraction 
comer  and  the  refleeting  surfaee. 

In  phase  two,  the  position  of  local  explosions  observed  on  the  refleeting  surface 
determine  the  range  of  aeeeptable  angles  and  offsets  for  reinitiation.  In  phase  three,  the  flame 
separation  and  shock  speed  after  the  second  diffraction  comer  determine  the  need  for  additional 
refleeting  surfaces  to  repeat  the  process  of  deeoupling  and  reinitiation  until  the  channel  height  is 
greater  than  the  critieal  channel  height. 

In  the  first  phase,  the  experimental  objectives  include  development  of  maps  of  the  flame 
separation  distance  and  shoek  strength  downstream  of  the  initial  diffraetion.  The  shape  of  the 
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first  diffraction  comer  dictates  the  rate  of  decrease  in  shoek  Maeh  number  and  inerease  in  shock- 
flame  separation.  To  investigate  the  effect  of  diffraction  angle  and  eomer  radius  on  shoek  deeay, 
a  seleetion  of  diffraction  angles  and  comer  radii  were  evaluated  (Fig.  5).  This  experimental 
sequenee  was  titled  “D”  as  a  shorthand  for  diffraetion.  The  experimental  results  inelude  eontour 
maps  of  shoek-flame  separation  distance  and  shoek  Mach  number  as  functions  of  diffraction 
angle,  comer  radius. 


Figure  5.  Initial  diffraction  design  parameters 

It  will  be  shown  that  a  eombination  of  high  shock  Mach  number  and  large  shock-flame 
separation  distance  are  preferred  for  the  highest  probability  of  a  local  explosion  without  trapping 
the  secondary  detonation. 

In  the  second  phase,  the  objeetive  is  to  obtain  the  loeation  of  loeal  explosions  on  the 
converging  ramp.  A  fixed,  first  diffraction  comer  geometry  will  keep  the  shoek  Mach  number 
and  shoek-flame  separation  profiles  eonstant  as  the  refleeting  surface  parameters  vary.  By 
manipulating  refleeting  surfaee  angle,  and  position  (horizontal  and  vertieal  distanee  to  the 
diffraetion  eomer)  one  will  systematieally  vary  the  strength  and  turning  angle  of  the  incident 
shoek  (Fig.  6).  The  goal  of  the  second  phase  was  to  find  an  optimum  location  and  angle  of  the 
refleeting  surface.  In  Ch.  Ill,  this  test  sequenee  is  labeled  “R”  as  a  shorthand  for  refleetion.  The 
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occurrence  and  position  of  loeal  explosions  define  the  design  spaee  for  the  reflecting  surface. 
Due  to  the  natural  statistieal  variation  in  the  diffracting  detonation,  the  shoek-flame  separation 
distanee  and  shoek  Mach  number  will  be  shown  to  have  large,  statistieal  variation.  As  a  result, 
the  oeeurrence  of  loeal  explosions  requires  a  statistieal  treatment.  The  results  in  Chapter  IV 
include  the  probability  of  loeal  explosion  as  a  funetion  of  reflecting  surface  angle,  vertieal 
distanee  from  the  diffraction  comer,  and  horizontal  distance  from  the  diffraetion  comer. 


Horizontal  Distanee,  xo  Loeal 


Figure  6.  Reflecting  surface  design  parameters 

The  purpose  of  the  third  phase  was  to  eomplete  the  transition  from  the  sub  eritieal 
channel  to  a  super  eritieal  channel  with  minimal  decoupling.  In  the  third  phase,  the  objeetive 
was  a  qualitative  examination  of  several  multi-reflection  geometries.  This  test  series  was  labeled 
“M”  in  CH.  Ill  as  shorthand  for  multiple  obstacles.  Deeoupling  after  the  first  reinitiation  of 
detonation  is  to  be  avoided  (Fig  7)  because  it  defeats  the  purpose  of  the  diffuser;  however  is  was 
universally  observed  in  the  R-series  test  eases. 

The  evolution  of  the  M-series  test  eases  was  based  on  the  observed  deeoupling  the 
previous  ease  beginning  with  the  most  consistently  reinitiating  case  from  the  R-series  (case  R2). 
In  case  Ml  the  obstacle  height  was  shorter  and  the  diffraction  comer  at  the  end  of  the  surfaee  had 
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a  larger  radius.  Later  the  number  of  reflecting  surfaces  was  increased  (M2-M5),  reflecting 
surfaces  were  added  to  the  initial  diffraction  corner  (cases  M7-M10),  and  the  decoupled  shock 
encountered  multiple  reflecting  surfaces  arranged  radially  (case  Mil). 

The  most  successful  geometry  in  this  research  utilized  multiple  reflection  surfaces 
interacting  separately  with  the  initially  decoupled  shock  (see  Fig.  33).  The  separately  reinitiated 
detonations  produced  local  explosions  in  the  region  downstream  of  the  reflecting  surfaces  due  to 
the  collision  of  diffracted  shockwaves  from  the  separate,  secondary  detonations  (see  Fig.  78). 
Spontaneous  local  explosions  are  the  defining  characteristic  of  detonation  diffraction  from  a 
critical  channel  (Soloukhin  and  Ragland,  1969).  The  occurrence  of  spontaneous  local  explosion 
in  a  case  where  the  initial  channel  height  is  subcritical  shows  an  improvement  gained  from  the 
addition  of  reflecting  surfaces. 


Figure  7.  Second  diffracting  corner  design  parameters 
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Units 


The  detonation  for  propulsion  eommunity  works  in  both  the  English  and  SI  unit  systems. 
To  appeal  to  a  broader  audienee  and  reduee  the  elutter  of  reporting  values  in  two  unit  systems, 
this  work  reports  only  the  SI  units. 

Organization 

This  dissertation  begins  with  a  detailed  examination  of  detonation  diffraction,  shock 
initiation,  and  detonation  kernel  development.  The  background  chapter  draws  from  relevant 
literature  to  describe  the  significant  phenomena  exploited  to  develop  a  detonation  diffuser.  The 
experimental  methodology  chapter  describes  experimental  methods,  measurement  techniques, 
equipment  requirements,  and  data  acquisition  systems.  The  analytical  methods  chapter  describes 
the  manual  and  automated  data  reduction  algorithms  and  the  associated  uncertainty.  The  results 
chapter  reports  the  observations  from  each  of  the  test  cases  and  describes  the  limits  of  the  design 
parameters.  A  conclusions  chapter  gathers  the  wisdom  gained  from  the  results  to  recommend  a 
functional  diffuser  geometry  and  additional  steps  toward  an  optimized  design.  Finally,  the 
bibliography  lists  the  literary  sources  used  throughout  the  paper. 
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II.  Background  and  Theory 


Overview 

This  chapter  details  the  relevant  portions  of  detonation  theory  and  empirical  evidence 
necessary  to  understand  the  reasoning  behind  the  experimental  methods  used  in  this  research. 
The  first  section  reviews  the  literature  concerning  diffraction  of  a  subcritical  detonation.  The 
second  section  examines  prior  work  on  both  normal  and  oblique  shock  reflections  leading  to 
detonation  initiation.  The  final  section  looks  at  the  development  of  a  detonation  kernel. 
Subcritical  Detonation  Diffraction 

Skews  (1967)  constructed  a  geometric  model  (Fig.  8)  for  the  head  of  a  disturbance 
propagating  into  the  fluid  behind  a  normal  shock  wave  during  diffraction.  The  Skews 
construction  is  useful  for  modeling  the  propagation  of  the  shockwave  as  it  decouples  from  the 
combustion  front.  It  lacks  any  treatment  of  heat  release  from  combustion  and  is  used  only  to 
model  the  shock  propagation  after  decoupling.  Figure  8  shows  the  state  of  the  diffracting  shock 
at  a  time,  At,  after  the  normal  shock  encounters  the  comer. 


Figure  8.  Construction  of  a  diffracting  shock  wave  (Skews  1967) 
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When  the  normal  shock  encounters  the  diffraction  comer,  an  expansion  wave  begins  at  the 
comer  and  traverses  the  normal  shock  at  the  post-shock  speed  of  sound  (a).  Meanwhile,  the 
unaffected  portion  of  the  shock  continues  to  propagate  at  the  original  velocity  (D).  The 
intersection  of  the  expansion  wave  and  the  shock  traces  a  straight  line  out  from  the  comer  at  an 
angle  a.  The  angle  depends  on  D,  a,  and  the  bulk  velocity,  u,  induced  by  the  shock  (Skews, 
1967). 

tan(a)  =  -  = - - -  (2) 

The  portion  of  the  shock  disturbed  by  the  expansion  wave  curves,  and  the  shock  velocity 
decreases.  When  the  initial  channel  height  is  subcritical  (see  Fig.  1),  the  reduced  shock 
compression  causes  the  shock  and  combustion  fronts  to  decouple  into  a  leading  shock  and  a 
trailing  flame.  When  the  charmel  height  is  greater  than  the  critical  height,  the  loss  of 
compression  is  insufficient  to  cause  global  decoupling  of  the  detonation. 

The  critical  channel  height  is  a  function  of  the  cross-section  of  the  channel,  and  the 
stability  of  the  detonation  wave  (Lee  2007).  Lee  (2007)  deemed  as  “unstable”  any  detonable 
mixture  that  resulted  in  a  detonation  with  cellular  stmcture.  For  the  2D  entrance  channel  used  in 
this  work,  Lee  found  that  the  critical  channel  height  was  six  times  the  cell  size  defined  as  A,. 

Pintgen  (2004)  examined  the  decoupled  shock  speed  and  shock-flame  separation  distance 
for  reactant  mixtures  of  hydrogen/oxygen/argon  and  hydrogen/N20  in  detonations  diffracting 
from  subcritical  channels.  Reported  were  the  shock  velocity  and  shock-flame  separation 
distance  as  functions  of  time  and  angle  measured  from  the  exit  plane  of  the  initial  channel, 
(labeled  |3  in  Fig.  9),  but  not  as  a  function  of  position.  Pintgen  used  four  methods  to  calculate 
distance  traveled  by  a  shock  between  photographic  frames.  The  first  used  the  measurement  from 
a  point  on  one  shock  to  the  closest  point  on  the  shock  in  the  following  frame,  the  "forward 
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closest"  method  (see  Fig.  9).  The  second  was  a  measurement  from  an  arbitrary  point  to  the 
closest  point  on  the  shock  in  the  preceding  frame  in  a  "backward  closest"  method.  The  third  and 
fourth  were  measurements  along  a  vector  normal  to  the  shock  front.  Measuring  to  the  following 
frame  was  "forward  normal"  and  measuring  to  the  preceding  frame  was  "backward  normal." 


Figure  9.  Measurement  methods  used  to  obtain  lead  shock  velocity  (Pintgen,  2004) 

The  velocity  profiles  obtained  from  these  measurements  indicate  where  the  shock  Mach  number 

decreases  most  slowly  giving  the  best  probability  of  reinitiation  when  the  shock  reflects. 

Pintgen  found  that  the  shock  speed  was  highest  near  the  centerline  of  the  channel  and 
decreased  as  p  increased  (Fig.  10).  The  shock  speed  also  decreased  rapidly  in  time  dropping 
over  50%  in  the  first  50  ps  of  diffraction.  Unfortunately,  Fig.  10  is  unsuited  for  the  purposes  of 
this  research  because  Pintgen  (2004)  considered  only  one  combination  of  diffraction  angle  and 
comer  radius  (90°  diffraction  angle,  0  mm  radius).  This  research  uses  a  configuration  with  the 
same  diffraction  angle  and  slightly  larger  diffraction  angle  (2.0  mm)  extensively  as  a  baseline  for 
comparison  to  different  diffraction  angle  and  comer  radii. 
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Angle  from  exit  plane  center  to  tube  ax  is  fdegl  Angle  from  exit  plane  center  to  tube  axis  fdeg] 


a)  2H2+02-h7Ar,  averaged. 


b)  H2-hN20,  averaged. 


c)  2H2+02-h7Ar.  raw  data.  d)  H2+N2O.  raw  data. 


Figure  10.  a)  Averaged  velocity  profiles  assuming  axis  symmetry  for  forward  and  backward  closest  point 
method,  2H2+02+7Ar,  Pq  =  1  bar.  Legend  gives  point  in  time  after  detonation  exited  the  tube,  b)  Averaged 
velocity  profiles  assuming  axis  symmetry  for  forward  and  backward  closest  point  method,  H2+N20  Pq  =  0.4 
bar.  Legend  gives  point  in  time  after  detonation  exited  the  tube,  c)  Normalized  velocity  obtained  with 
forward  closest  point  technique  for  2H2+02+7Ar,  Pq  =  1  bar.  d)  Normalized  velocity  obtained  with  forward 
closest  point  technique  for  H2+N20  Pq  =  0.4  bar  (Pintgen,  2004). 


A  lower  diffraction  angle  alone  can  prevent  decoupling  (Nettleton  1987).  In  Fig.  10,  the 
diffraction  angle  was  90"^.  Nettleton  (1987)  predicted  a  maximum  diffraction  angle  below  which 
the  rate  of  expansion  is  small  enough  to  prevent  decoupling  of  the  detonation.  For  a  reactant 
mixture  with  y  =  1.4,  Nettleton  (1987)  predicts  no  decoupling  at  diffraction  angles  below  14. 5"^. 
Nettleton  (1987)  derived  the  minimum  angle  from  a  combination  of  the  Chester-Chisnell- 
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Whitham  shock  diffraction  theory  and  the  Chapman-Jouguet  detonation  wave-speed  theory. 
According  to  Nettleton’s  analysis,  the  minimum  angle  maintains  a  sufficient  shock  velocity 
along  the  diverging  wall  for  detonation.  Nettleton  (1987)  did  not  incorporate  triple  point 
interactions  that  aid  in  sustaining  detonations,  and  experimental  work  by  the  author  shows  that 
some  decoupling  occurs  for  angles  as  small  as  14°  (Stevens  at  al.,  2011(a))  and  at  15°  a 
detonation  does  not  fully  decouple  (see  Fig.  53). 

The  impetus  to  vary  the  diffraction  comer  radius  in  the  current  work  was  prompted  by 
observations  in  crossover  tubes  by  Nielsen  et  al.  (2011)  who  varied  the  crossover  tube  geometry 
and  found  a  delay  in  decoupling  when  the  diffraction  comer  radius  was  increased.  A  25.4  mm  in 
comparison  to  a  2.0  mm  radius  indicated  a  qualitative  delay  in  decoupling  of  the  incident 
detonation  front.  Nakayama  et  al.  (2012)  also  reported  increased  wave  speed  as  the  inner  radius 
increased  in  a  curved,  square  cross-section  channel.  In  Fig.  11,  Nakayama  et  al.  show  the  onset 
and  increase  of  decoupling  as  the  inner  wall  radius  decreases.  The  minimum  ratio  of  radius  to 
cell  size  was  21  for  prevention  of  the  unstable  mode  where  decoupling  occurs.  For  a 
stoichiometric  mixture  of  hydrogen  and  air,  the  cell  size  is  8.19  mm  (Ciccarelli  et  al.,  1994),  and 
thus  the  minimum  radius  to  prevent  decoupling  according  to  Nakayama  et  al.  would  be  172.2 
mm.  The  current  work  differs  because  there  is  no  outer  wall  to  reflect  triple  points  and  restore 
detonation  in  the  unstable  mode.  It  is  unknown  what  the  minimum  radius  would  be  without  the 
outer  wall,  and  Nakayama  is  the  best  available  prediction.  Because  decoupling  is  desirable  to 
separate  the  shock  and  flame  prior  to  reinitiation,  the  diffraction  comer  radius  was  always  less 
than  the  predicted  172.2  mm  in  the  current  work.  Two  diffraction  comer  radii  were  studied  to 
determine  the  sensitivity  of  decoupling  to  the  radius.  The  next  section  explains  why  decoupling 
is  necessary  to  make  the  transition  from  a  subcritical  channel  to  a  supercritical  channel. 
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Figure  11.  Diffraction  in  channels  of  decreasing  inner  radius.  Stoichiometric  mixture  of  ethylene  and 
oxygen.  Decoupling  is  visible  when  ri  <  40  mm  in  panels  c,  d,  and  e.  (Nakayama  et  al.,  2012) 
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Detonation  Initiation  via  Shock  Reflection 

Localized  explosions  are  the  preeursor  to  stable,  cellular  detonation.  Local  explosions 
occur  in  both  DDT  (Urtiew  and  Oppenheim,  1966)  and  diffracting  detonations  exiting  critieal 
height  channels  (Soloukhin  and  Ragland,  1965).  In  DDT,  loeal  explosions  occur  in  the  spaee 
between  the  leading  normal  shock  and  the  aeeelerating  combustion  front  (Urtiew  and 
Oppenheim,  1966).  In  critieally  diffraeting  detonations,  loeal  explosions  oeeur  near  the  in  the 
spaee  between  the  decoupled  shoek  and  eombustion  front  (Soloukhin  and  Ragland,  1965).  In 
both  situations,  the  loeal  explosion  results  in  a  detonation  kernel  that  grows,  develops  cellular 
structure,  and  stabilizes  in  a  channel  to  become  a  planar  detonation  front. 

Localized  explosions  also  oeeur  when  high  Mach  number  shoeks  refleet  from  surfaees 
(Brown  and  Thomas,  2000  and  Thomas  et  al.,  2002).  Early  detonation  experiments  found  that 
detonation  eould  also  be  initiated  by  normal  shock  reflection  at  the  end  wall  of  a  shoek  tube 
(Saitsev  and  Soloukhin,  1958  and  Strehlow  and  Cohen,  1962).  In  either  case,  the  formation  of  a 
detonation  depends  on  local  speed  of  sound,  reflecting  surface  area,  and  induetion  delay  of  the 
detonable  mixture  (Thomas  et  al.,  2002).  A  critieal  eondition  below  which  detonation  initiation 
is  improbable  is  p  <  1,  where  r\  is  the  ratio  of  surface  height  (h)  to  the  produet  of  post-shoek 
speed  of  sound  (a)  and  induction  delay  (x)  shown  in  Eq.  3. 

—  =  V  (3) 

Induetion  delay  is  the  time  that  passes  between  the  shoek  refleetion  and  the  onset  of  heat  release 
by  chemical  reaction  of  the  reactants.  Thomas  et  al.  defined  surfaee  height  for  a  normal  shock 
reflection.  In  the  eurrent  work  the  definition  has  been  generalized  for  oblique  reflections  as  the 
perpendieular  distanee  from  the  channel  wall  to  the  end  of  the  obstruction  (Fig.  12.). 


17 


Reflecting  Surfaces 


Figure  12.  Surface  height  for  normal  and  oblique  reflecting  snrfaces 


To  address  the  suggestion  by  Brown  and  Thomas  that  a  shock/boundary  layer  interaction 
might  be  necessary  to  initiate  detonation,  Thomas  et  al.  (2002)  utilized  a  normal  shock  reflection 
from  the  end  of  a  cylinder  to  quantify  the  critical  conditions  (Fig.  13).  Prior  to  the  shown 
frames,  a  right  traveling  shock  impacted  the  circular  face  of  the  cylinder  causing  a  reflected 
shock  and  a  flame  front  to  form.  For  the  cylinder,  Thomas  et  al.  substituted  the  cylinder  radius 
for  the  surface  height  to  adapt  Eq.  to  the  new  configuration.  In  the  experimental  frames  on  the 
left,  q  was  0.93  and  the  mixture  ignited  after  the  reflection,  but  did  not  detonate.  In  the 
simulated  frames  on  the  right,  q  was  1.00  and  the  shock  reflection  initiated  detonation.  Taken 
together,  the  experimental  and  simulation  results  validate  the  substitution  of  radius  for  height  in 
the  critical  condition  for  detonation  initiation  and  remove  the  boundary  layer  requirement 
suggested  by  Brown  and  Thomas.  The  current  work  relies  on  this  finding  first  when  extending 
the  definition  of  surface  height  to  oblique  reflecting  surfaces  and  later  when  using  multiple 
reflecting  surfaces  to  initiate  multiple,  separate  detonations  from  the  same  decoupled  shock. 
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Figure  13.  Detonation  initiation  by  reflecting  normal  shock  (Thomas  et  al.,  2002) 


The  earliest  evidence  of  reinitiation  by  oblique  shock  reflection  came  from  crossover  tube 


studies  (Nielsen  et  ah,  2011).  A  crossover  tube  also  causes  diffraction  of  the  detonation,  and 


reinitiation  can  occur  when  the  decoupled  shock  reflects  from  the  reflecting  surface  (Fig.  14). 
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Figure  14.  Shock  initiation  in  a  crossover  tube  (Nielsen  et  al.,  2011) 

The  decoupled  shock  was  nearly  normal  to  the  reflecting  surface  on  impact  and  the  secondary 

detonation  formed  in  the  crossover  tube  did  not  propagate  back  to  the  primary  tube  because  the 
combustion  front  consumed  the  reactants  before  the  detonation  arrived.  Diffraction  at  the  exit  of 
the  crossover  tube  caused  the  secondary  detonation  to  decouple  and  it  was  necessary  to  add  DDT 
obstacles  to  the  secondary  tube  to  transition  to  detonation  (Nielsen,  2011). 

The  crossover  geometry  in  Fig.  14  has  an  identical  diffraction  comer,  but  the  reflection 
surface  is  perpendicular  to  the  channel  wall  and  not  attached.  The  result  is  a  secondary 
detonation  that  is  isolated  or  from  the  remaining  reactants  by  combustion  products  or  “trapped” 
(Fig.  16).  The  diffraction  angle  is  90°  and  the  corner  radius  is  2.0  mm.  The  reflection  angle  is 
90°  and  the  leading  edge  of  the  reflecting  surface  is  38.1  mm  downstream  of  the  diffraction 
comer.  The  vertical  offset  is  0  mm.  Local  explosions  occur  after  the  decoupled  flame  front 


20 


reaches  the  reflecting  surface  preventing  the  secondary  detonation  from  propagating 
downstream.  Trapped  detonations  are  undesirable  in  a  detonation  diffuser  because  a  trapped 
detonation  cannot  propagate  downstream.  The  key  to  preventing  a  trapped  secondary  detonation 
is  minimizing  the  delay  between  shock  reflection  and  reinitiation  known  as  the  induction  delay 
(x)  while  maximizing  the  separation  distance  between  the  decoupled  shock  and  flame  to  allow 
the  secondary  detonation  time  to  propagate  downstream  before  the  decouple  combustion  front 
arrives. 

The  crossover  tubes  used  in  detonation  branching  configurations  (Nielsen,  2011  and 
Camardo,  2012)  share  phenomena  with  the  detonation  diffuser.  The  decoupling  and  reinitiation 
processes  are  identical.  Stable  detonation  waves  undergo  diffraction  at  the  crossover  and 
decouple  (Fig.  15b).  Shock  reflection  on  the  opposite  wall  of  the  crossover  tube  reinitiates 
detonation  (Fig.  15c),  and  diffraction  occurs  again  at  the  end  of  the  crossover  tube  (Fig.  15d). 
Because  of  the  similarities,  crossover  tube  studies  were  a  good  starting  point  for  the  current 
diffuser. 
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a)  Frame  0:  0  |as 
Second  diffraction  com 
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First  diffraction  comer 
Planar  primary  detonation 

b)  Frame  1:  38.0  ps 
Decoupling 
Expansion  head 
Unaffected  detonation 


c)  Frame  3:  76  ps 
Secondary  detonation 
Decoupled  shock 
Decoupled  combustion 


d)  Frame  4:  144  ps 
Decoupled  shock 
Decoupled  combustion 


Figure  15.  Diffraction  and  reinitiation  in  a  crossover  tube  (Nielsen  et  al.,  2011) 
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The  most  significant  difference  between  a  crossover  tube  and  a  diffuser  is  the  orientation 
of  the  reflecting  surface.  The  reflecting  surface  is  typically  parallel  to  the  initial  detonation  front 
in  a  crossover  tube  to  reduce  the  volume  of  the  crossover  tube.  Flow  loss  considerations  are 
secondary  to  minimizing  the  volume  of  fuel/air  mixture  needed  to  fill  the  crossover,  as  the 
crossover  tube  produces  little  thrust.  In  a  detonation  diffuser,  the  diffuser  volume  contributes  to 
thrust  and  the  flow  losses  are  more  prevalent  than  in  a  crossover  configuration.  As  a  result,  the 
desired  reflection  angle  is  as  small  as  reliable  reinitiation  allows.  The  crossover  tube  studies 
suggested  that  increasing  the  radius  of  the  initial  diffraction  comer  delays  decoupling  within  the 
crossover  tube  (Fig.  16).  In  the  second  frame  of  Fig.  16a,  the  separation  distance  is  4.89  ±  0.98 
mm,  and  in  Fig.  16b,  the  separation  distance  is  4.03  ±  0.98  mm.  The  difference  was  small  and 
the  uncertainties  were  large  enough  that  a  better  experiment  was  needed  to  draw  a  statistically 
significant  conclusion,  but  the  trend  is  encouraging  so  diffraction  comer  radius  was  included  as  a 
parameter  in  the  D-series  test  cases.  Unlike  the  separation  distance  difference,  shock  reflection 
induced  detonations  were  obvious  in  the  crossover  tube  videos. 


Figure  16.  Delayed  decoupling  due  to  large  corner  radius  (Nielsen  et  al.,  2011) 
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Reinitiation  was  observed  on  both  flat  and  concave  reflecting  surfaces  (Fig.  17).  The 
increased  diffraction  angle  due  to  the  concavity  ensured  decoupling  in  the  second  detonation 
tube.  The  remaining  test  cases  avoided  concave  reflecting  surfaces  to  reduce  diffraction  at  the 
end  of  the  reflecting  surfaces. 


a)  Flat  wall  b)  25  mm  concave  radius 


Figure  17.  Local  explosion  due  to  shock  reflection  from  flat  and  concave  surfaces  (Nielsen  et  al.  2011) 
Oblique  reflections  also  initiate  detonation  provided  the  compression  is  sufficient. 

Detonation  initiation  after  an  unsteady  oblique  shock  reflection  was  observed  by  Stevens  et  al. 
(2011).  In  Fig.  18,  Stevens  et  al.  compare  a  2D  simulation  of  diffraction  and  reinitiation  to 
experimental  frames  from  a  38  cm  high,  50  cm  deep  channel.  In  the  first  frame  of  the 
simulation,  diffraction  begins,  but  it  is  not  yet  visible  in  the  experiment.  In  the  second  frame  of 
both  series,  separation  is  visible.  In  the  simulation,  the  separation  region  is  a  light  blue  region 
between  the  shock  and  the  flame.  In  the  experiment,  the  separation  region  is  a  subtle  dark  band 
following  the  shock  front.  In  the  last  frame  of  each  series,  detonation  has  reinitiated  along  the 
reflecting  surface. 
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Figure  18.  Detonation  initiated  by  obliqne  reflection  (Stevens  et  al.  2011) 


The  geometry  in  Fig.  18  contains  most  but  not  all  of  the  parts  of  a  detonation  diffuser.  The 
configuration  has  a  diffraction  angle  of  90°,  a  comer  radius  of  0.0  in  the  simulation  (2.0  mm  in 
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the  experiment),  a  reflection  angle  of  14°  and  a  vertical  offset  of  -12.7  mm.  There  is  no  visible 
diffraction  at  the  end  of  the  reflection  surface  because  the  end  is  not  visible. 

Detonation  Kernel  Development 

In  Figure  15,  the  local  explosion  occurs  near  the  midpoint  of  the  reflecting  surface.  Since 
there  is  significant  shot-to-shot  variation  expected  in  the  decoupled  shockwave,  an  excess  of 
reflecting  surface  will  be  necessary  to  improve  the  probability  of  a  local  explosion,  thus  the 
detonations  begun  by  local  explosions  in  a  detonation  diffuser  need  to  propagate  along  the 
remainder  of  the  reflecting  surface  and  past  the  second  diffraction  comer  at  the  end  for  the 
diffuser  to  be  effective.  The  natural  evolution  the  local  explosion  into  a  detonation  wave 
determines  the  speed  and  shape  of  the  detonation  front  as  it  encounters  the  diffraction  comer  at 
the  end  of  the  surface.  The  secondary  detonation  will  initially  be  overdriven  (Schauer  et  al., 
2005),  and  the  excess  wave-speed  will  be  important  to  overcome  the  loss  due  to  diffraction  at  the 
second  diffraction  comer.  To  determine  the  maximum  angle  of  diffraction  that  the  detonation 
can  tolerate  without  decoupling,  it  is  necessary  to  characterize  the  evolution  of  the  detonation 
from  local  explosion  to  the  quasi-stable  cellular  mode  exhibited  by  detonations  in  channels 
(Urtiew  and  Oppenheim  1966).  In  the  cellular  mode,  the  local  detonation  Mach  number  varies  in 
a  repeating  pattern  dependent  on  the  propagation  of  transverse  shock  waves  that  intersect  the 
detonation  front.  The  cellular  mode  is  quasi-stable  because  the  local  wave  speed  varies,  but  the 
average  speed  of  advance  does  not. 

Experimental  studies  of  DDT  show  initial  detonation  kernel  propagation  dominated  by 
the  local  speed  of  sound  and  shock  induced  velocity.  Urtiew  and  Oppenheim  (1966)  captured 
the  propagation  in  a  25  mm  x  38  mm  rectangular  cross-section  channel  with  framing  schlieren 
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images  (Fig.  19).  In  the  first  frame  of  Fig.  19,  the  local  explosion  has  occurred  in  the  boundary 
layer  along  the  top  wall  of  the  channel  sometime  between  30  and  35  ps  (between  consecutive 
frames).  The  detonation  convected  with  the  bulk  motion  behind  the  normal  shock,  and  expanded 
until  the  detonation  front  caught  up  with  the  normal  shock  at  approximately  the  40  ps  mark.  Part 
of  the  wave  front  continued  to  propagate  through  the  post  shock  region  and  encountered  the 
deflagration  front  at  the  same  time.  The  detonation  front  reached  the  bottom  wall  of  the  channel 
consuming  the  reactants  before  the  deflagration  arrived  just  after  the  45  ps  mark.  Without 
reactants,  the  deflagration  perished  by  50  ps.  The  part  of  the  detonation  that  passed  the  normal 
shock  slowed  because  of  the  lower  pressure  and  temperature  ahead  of  the  shock.  The  detonation 
front  picked  up  velocity  instabilities  from  interaction  with  normal  shock.  The  instabilities  grow 
in  size  from  40  ps  to  50  ps  as  the  detonation  front  continues  downstream  and  will  eventually 
stabilize  as  cellular  structure  (Urtiew  and  Oppenheim,  1966). 
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Figure  19.  DDT  observed  in  hydrogen/oxygen  (Urtiew  and  Oppenheim,  1966). 
Structures  highlighted  by  the  author 


Local  explosions  produce  a  smooth  spherical  blast  wave  that  evolves  cellular  structure 
after  developing  instability  (Gamezo  et  ah,  1999).  Numerical  simulation  of  blast  waves  (Fig.  16) 
expands  the  time  scale  of  the  onset  of  instability  showing  the  transition  from  a  smooth  blast  wave 
to  quasi-stable  cellular  detonation  that  happens  too  quickly  to  capture  with  the  imaging  technique 
used  by  Urtiew  and  Oppenheim  in  Fig.  19. 
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Figure  20.  Evolution  of  cellular  structure  from  a  smooth  blast  wave  (Gamezo,  1999) 

Gamezo  et  al.  (1999)  found  that  triple  points  form  as  the  result  of  instability  in  the  blast 

wave  and  weak  turbulence  in  the  reactants.  The  perturbations  cause  wrinkling  in  the  smooth 
wave,  and  transverse  waves  form.  At  L  =  6  mm  in  Fig.  20,  the  cell  size  is  small  because  there 
are  many  perturbations  of  the  blast  wave.  As  the  detonation  wave  slows,  some  triple  points 
disappear  and  the  cell  size  increases.  Over  time,  the  cellular  detonation  wave  slows 
asymptotically  to  the  CJ  speed  and  the  cell  size  becomes  constant. 

In  the  detonation  diffuser,  secondary  detonations  that  form  along  the  reflecting  surface  do 
not  encounter  the  far  wall  of  the  channel  before  the  end  of  the  reflecting  surface,  and  the 
collected  data  will  show  whether  the  secondary  detonation  or  the  decoupled  shock  will  reach  the 
end  of  the  second  diffraction  comer  first.  How  the  secondary  detonation  reacts  to  the  diffraction 
comer  at  the  end  of  the  surface  will  depend  on  the  shock  speed  and  number  of  triple  points 
present  at  that  moment.  Since  no  published  data  describes  that  interaction,  this  research  reports 
the  speed  of  secondary  detonations  before  and  after  encountering  the  second  diffraction  comer 
and  the  shock  Mach  number  and  shock-flame  separation  distance  when  the  secondary  detonation 
decoupled. 


29 


The  literature  coneeming  detonation  diffraction,  shock  initiation  of  detonation,  and 
evolution  of  a  detonation  front  from  local  explosions  suggests  three  sets  of  measurements 
necessary  to  enable  analysis  of  a  detonation  diffuser.  The  first  measurements  are  the  shock 
Mach  number  and  shock-flame  separation  distance  as  functions  of  diffraction  angle  (0i)  and 
comer  radius  (rj).  The  second  measurements  are  location  and  probability  of  local  explosions  as 
functions  of  reflecting  surface  angle  and  position.  The  final  set  of  measurements  are  shock  speed 
and  separation  distance  after  diffraction  of  a  secondary  detonation.  This  research  reports  all 
three  sets  experimentally  as  any  numerical  solution  would  require  experimental  validation.  The 
next  section  describes  the  experimental  methods  employed.  A  successful  detonation  diffuser 
initiates  one  or  more  secondary  detonations  that  do  not  decouple  when  the  final  channel  height  is 
greater  than  the  critical  channel  height. 
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III.  Experimental  Methodology 

Overview 

The  experimental  methods  in  this  research  fall  under  one  of  two  headings:  experimental 
techniques  or  data  collection  methods.  Experimental  techniques  describe  the  equipment  and 
procedures  for  initiating  detonation  and  subjecting  detonations  to  test  articles.  Data  collection 
methods  describe  the  apparatus  and  procedures  that  collect  information  from  the  detonation  as  it 
interacts  with  test  articles.  Chapter  IV  will  discuss  data  reduction  and  uncertainty.  The 
techniques  and  collection  methods  depend  on  each  other  and  on  the  characteristics  of  a 
detonation  wave.  As  a  result,  a  combination  of  careful  planning  and  mid-course  adjustments 
ensured  relevant,  accurate  measurements. 

Experimental  Techniques 

Investigating  detonation  diffraction  and  shock  initiation  requires  a  source  of  repeatable 
detonations.  The  research  PDE  at  the  Detonation  Engine  Research  Facility  is  one  such  source. 
Schauer  et  al.  (2001)  first  published  the  details  of  the  engine.  Nielsen  (2011)  included  a 
description  of  the  configuration  used  for  schlieren  visualization  of  crossover  geometries  (Fig. 
21).  The  PDE  head  and  Tube  2  were  kept  and  a  new  section  consisting  of  a  narrow  channel  with 
a  large,  instantaneous  increase  in  channel  height  (Fig.  22)  replaced  Nielson's  "test  rig"  and 
"manifold"  sections.  Hence,  an  expansion  section  was  introduced  in  the  upward  direction,  which 
is  physically  much  like  Fig.  4  except  that  the  expansion  is  in  the  opposite  direction. 
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Figure  21.  PDF  as  configured  for  crossover  study  (Nielsen  et  al.,  2011) 
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Figure  22.  CAD  model  of  optical  test  section 

The  PDE  operates  on  a  wide  variety  of  gaseous  and  liquid  fuels  offering  a  wide  range  of 
cell  sizes.  The  detonation  frequency  is  adjustable  from  8  Hz  to  40  Hz,  and  the  ignition  can  be  set 
to  operate  in  “burst  mode”  firing  for  a  predetermined  number  of  cycles.  Adjusting  the 
equivalence  ratio  of  the  mixture  gives  control  of  cell  size  for  a  specific  fuel/air  mixture 
(Ciccarelli  et  ah,  2004). 

The  design  of  the  research  PDE  presents  some  undesired  effects.  The  detonation  tube 
pressure  at  ignition  is  not  explicitly  controlled,  and  the  dynamics  of  filling  cause  variations  of 
pressure  throughout  the  detonation  tube  (Helffich,  2006).  Applying  an  ignition  delay  of  4  ms 
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after  the  elose  of  the  fill  valves  mitigated  the  variations.  The  periodic  fill  process  also 
contributes  to  local  variations  in  equivalence  ratio.  Fuel  flows  into  the  airstream  constantly  as 
the  air  flow  varies  resulting  in  locally  rich  and  lean  conditions.  Use  of  a  second  detonation  tube 
180°  out  of  phase  with  the  test  section  (Tube  4  in  Fig.  21)  halved  the  time  that  the  fill  valves 
were  closed  and  reduced  the  variations  so  that  they  were  less  than  5%.  With  the  mitigation  of  its 
undesired  characteristics,  the  research  PDF  served  as  the  source  of  reactants  and  detonations  for 
all  test  cases. 

After  fill,  ignition,  and  DDT,  detonations  passed  into  an  optical  test  section  containing 
the  various  test  articles  (Fig.  22).  An  adapter  gradually  transitioned  from  the  circular  detonation 
tube  cross-section  to  the  narrow  channel  in  the  test  section.  The  channel  was  6.35  mm  wide 
(from  window  to  window)  or  77%  of  the  8.19  mm  cell  width  for  stoichiometric  H2/Air  at 
atmospheric  pressure  and  temperature  (Ciccarelli  et  ah,  1994).  The  small  width  of  the  channel 
suppressed  cellular  structure  in  that  dimension.  The  optical  section  begins  with  a  section  of 
subcritical,  rectangular  channel  (h/A,  =  6.2)  and  opens  into  a  taller  section  containing  the 
geometry  under  study  (Fig.  22).  A  stoichiometric  mixture  of  H2/Air  was  used  in  all  test  cases 
for  a  consistent  ratio  of  initial  channel  height  to  cell  size. 

After  the  adapter,  the  channel  size  was  constant  for  127  mm  allowing  the  wave  speed  to 
stabilize.  The  entrance  channel  height  was  50.8  mm  opening  up  to  a  maximum  height  of  191 
mm.  The  test  section  had  optical  access  for  schlieren  visualization  via  two  polycarbonate 
windows.  The  windows  were  each  12.7  mm  thick  and  tolerated  the  impulsive  detonation  loading 
without  any  evidence  of  fatigue  after  hundreds  of  detonations. 

An  unexpected  phenomenon  encountered  early  in  testing  was  the  propagation  of  strain 
waves  through  the  polycarbonate  windows.  The  strain  caused  a  small  change  in  the  refraction  of 
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light  through  the  window.  The  resulting  distortion  of  the  light  was  on  the  order  of  the  refraction 
caused  by  density  gradients  within  the  channel,  and  appeared  as  light  and  dark  bands  in  the 
recorded  images  (Fig.  23).  A  strain  gage,  attached  to  an  expended  window,  verified  that  the 
waves  were  due  to  the  bending  of  the  windows  and  not  due  to  density  gradients  within  the  test 
section.  Initially,  a  software  filter  attempted  to  remove  the  waves  from  the  recorded  images,  but 
proved  unable  to  discern  the  strain  waves  from  flames.  Reducing  the  sensitivity  of  the  schlieren 
system  reduced  the  visibility  of  the  strain  waves  with  acceptable  results,  and  the  affected  cases 
were  repeated. 


Figure  23.  Visible  strain  waves 

The  test  section  survived  the  intense  pressures  and  temperatures  associated  with 
detonation.  Peak  pressures  in  excess  of  3  MPa  and  peak  temperatures  near  3000  K  are  typical  of 
a  detonation  front  (Zeldovich,  1956),  but  the  maximum  values  are  short  lived,  and  12.7  mm  thick 
polycarbonate  was  an  acceptable  material  for  windows.  Brittle  surface  coatings  intended  to 
improve  scratch  resistance  were  tried,  but  the  strain  on  the  windows  caused  the  coating  to 
fracture.  Without  a  scratch-resistant  coating,  the  windows  regularly  suffered  abrasion  from  the 


34 


test  articles.  The  heat  from  local  explosions  did  bum  away  a  small  amount  of  the  window 
surface  after  several  cycles  (Fig.  24),  but  the  damage  from  abrasion  was  far  more  significant. 
The  windows  were  replaced  whenever  the  damage  obscured  the  shock  and  flame  within  the 
channel. 


Figure  24.  Scorch  marks  (red  arrows)  on  a  polycarbonate  window. 


The  test  section  allowed  for  simple  exchange  of  test  articles.  Quick  changeover  was 
important,  as  testing  time  is  limited  at  the  DERF.  Two  studs  permitted  the  removal  of  one 
window  at  a  time.  The  test  article  rested  on  the  remaining  window  until  secured.  Safety  wire 
secured  the  test  article  to  the  side  of  the  test  section  until  the  bolts  holding  the  windows  were 
tight,  preventing  movement  of  the  test  article. 

Data  Collection  Methods 

Schlieren  visualization  was  the  preferred  measurement  technique  for  detailed  study  of 
diffraction  and  local  explosions.  Schlieren  visualization  depicts  density  gradients  making  both 
shock  and  combustion  fronts  visible  (Fig.  25).  A  minimal  schlieren  visualization  system  is 
composed  of  a  light  source,  two  focusing  mirrors,  a  knife-edge,  and  a  screen  on  which  to  project 
the  image.  A  camera  usually  replaces  the  screen  to  capture  video  of  unsteady  flows.  Appendix 
A  describes  the  system  in  full.  Because  of  its  ability  to  visualize  both  shock  and  combustion  and 
to  record  at  high  frame  rates,  schlieren  visualization  saw  extensive  use  throughout  the  diffuser 
development. 
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Figure  25.  Schlieren  image  of  decoupling  detonation 

The  schlieren  technique  is  analog,  and  only  the  camera  limits  frame  rate  and  resolution. 

Cellular  structure  and  wave  speed  dictate  the  resolution  and  frame  rate  requirements.  The 
minimum  resolution  (measured  in  mm/pixel)  necessary  to  image  cellular  structure  is  half  the  cell 
width  per  pixel  according  to  the  Nyquist-Shannon  sampling  theorem,  and  the  maximum  temporal 
resolution  was  240,000  frames  per  second.  All  images  had  a  spatial  resolution  of  0.650 
mm/pixel.  The  temporal  restriction  relaxes  for  diffracting  detonations,  as  it  is  unnecessary  to 
capture  the  transverse  wave  collisions.  Instead,  the  temporal  resolution  depends  on  the 
dimensions  of  the  test  article.  Resolutions  in  this  work  ranged  from  4.76  ps  (210000  frame/s)  to 
25.0  ps  (40000  frame/s)  often  with  data  from  multiple  frame  rates  combined  into  one  data  set. 
The  exposure  of  each  image  was  293  ns,  and  the  resulting  motion  blur  at  Chapman  Jouguet 
speed  (-1971  m/s)  was  0.578  mm.  The  RMS  uncertainty  due  to  motion  and  spatial  resolution 
was  0.871  mm. 

Chemiluminescence  is  a  verification  tool  for  identifying  combustion  separately  from 
other  structures  in  the  flow.  A  camera  recorded  images  of  the  chemiluminescent  emissions  from 
chemical  reactions  for  a  crossover  tube  geometry  in  experiments  carried  out  in  prior  research 
(Fig.  26)  (Nielsen,  2011).  Because  the  same  camera  was  used  to  record  the  schlieren  and 
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chemiluminescence  images,  the  same  resolution  and  uncertainty  considerations  apply  to  both. 
No  equipment  beyond  a  zoom  lens,  a  high  frame  rate  camera,  and  optical  access  to  the 
detonation  are  needed,  making  setup  faster  than  schlieren  visualization.  The  weakness  of 
chemiluminescence  is  its  inability  to  detect  non-emitting  phenomena,  and  thus  the  decoupled 
shock  that  initiates  a  secondary  detonation  is  invisible.  Because  it  added  no  new  information, 
chemiluminescence  was  only  as  a  verification  tool  for  flame  from  propagation  and  local 
explosion  detection. 


Figure  26.  Schlieren  (left)  and  Chemiluminescence  (right)  of  decoupling  detonation  (Nielsen,  2011) 

Ion  probes  are  a  simple  solution  for  measuring  wave  speed  when  optical  access  is  not 
available.  Ion  probes  are  essentially  capacitors  that  close  a  circuit  when  combustion  ions  are 
present.  Wave  speed  is  computed  as  linearly  proportional  to  the  time  of  flight  between  two 
probes.  Two  pairs  of  probes  measured  wave  speed  upstream  of  the  optical  section  to  verify 
detonation  prior  to  entering  the  optical  test  section.  Low  difference  between  the  wave  speeds  at 
each  location  confirmed  that  the  detonation  was  propagating  at  the  CJ  speed.  The  uncertainty  in 
wave  speed  measured  by  ion  probes  is  a  function  of  the  distance  between  the  probes,  the 
sampling  frequency,  and  interpretation  of  the  capacitor-like  waveform.  The  sample  rate  for  ion 
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probes  in  this  work  was  a  constant  1  MHz,  and  the  distance  between  adjacent  ion  probes  was 
always  150  mm.  The  uncertainty  in  distance  between  the  probes  was  0.8  mm,  the  uncertainty  in 
time  is  0.5  ps,  and  the  uncertainty  in  the  waveform  was  12.8  m/s.  The  total  uncertainty  in  wave- 
speed  was  30.0  m/s. 

Calibration 

Shock  speed  and  shock-combustion  separation  distance  derive  from  the  position  of  the 
shock  and  flames  in  each  image.  To  measure  the  position  of  these  structures,  image  calibration 
was  necessary.  The  position  was  a  function  of  the  optical  magnification,  pixel  size,  and  distance 
from  the  camera  to  the  test  section.  To  bypass  deconvolving  the  effects  of  all  three  variables, 
pixel  positions  were  calibrated  using  the  initial  channel  height.  The  initial  channel  was  always 
visible,  and  it  had  a  known  height  of  50.8  mm.  The  calibrated  size  and  center-to-center  distance 
of  the  square  pixels  averaged  0.647  mm  (Fig.  27),  and  the  uncertainty  in  position  across  all  cases 
was  0.324  mm. 


0.647 


- > 


Figure  27.  Calibrated  pixel  dimensions 


Test  Cases 

This  research  set  out  to  accomplish  three  phases  of  study.  The  first  phase  will  quantify 
the  shock  Mach  number  and  separation  distance  after  a  detonation  encounters  a  diffraction 
comer  (D-series).  The  second  will  determine  the  reflecting  surface  angle  and  position  that  offer 
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the  best  chance  of  reinitiation  (R-series).  The  third  phase  will  investigate  multiple  obstacle 
geometries  to  narrow  the  field  of  possible  configurations  (M-series).  The  test  cases  mirror  the 
phases  of  the  research. 

There  were  three  sets  of  test  cases  (diffraction  “D”,  reflection  “R”,  and  multiple  obstacles 
“M”)  as  well  as  a  single  control  case.  This  section  contains  dimensioned  drawings  of  each.  The 
control  case  was  a  straight  channel  the  same  height  as  the  entrance  of  the  test  section  (50.8  mm) 
and  300  mm  in  length  (Fig.  28). 


Because  each  case  included  a  unique,  handmade  test  article,  the  test  matrix  was  kept  coarse 
while  still  bounding  the  limits  of  decoupling  or  reinitiation,  except  in  the  multiple  obstacle  cases. 
Those  latter  tests  included  a  series  of  incremental  changes  intended  to  approach  successful 
transition  to  a  super-critical  channel. 

In  the  diffraction  cases,  two  parameters  were  examined:  diffraction  angle  (0)  and  comer 
radius  (r)  (Fig.  28).  Table  1  lists  the  four  configurations  tested.  Case  D1  was  a  control  with  no 
diffraction.  Data  from  the  control  established  the  CJ  speed  for  the  remaining  cases.  Cases  D2 


39 


and  D3  compared  diffraction  angles  while  keeping  eomer  radius  constant.  Decoupling  was 
expected  and  oeeurred  in  Case  D2,  but  in  Case  D3  the  diffraetion  angle  was  0.5°  greater  than  the 
theoretieal  limit  proposed  by  Nettleton  (1983).  The  presence  of  transverse  shockwaves  not 
included  in  the  Nettleton’ s  analysis  made  it  unelear  if  decoupling  would  occur.  The  result  was 
partial  deeoupling  which  will  be  diseussed  in  Chapter  IV.  Cases  D2  and  D4  varied  the  diffraetion 
comer  radius  while  keeping  the  angle  constant,  hr  Case  D2,  the  radius  was  small  enough  for  the 
expansion  to  be  considered  instantaneous.  In  Case  D4,  the  radius  was  the  same  as  in  the  Nielsen 
et  al.  (201 1)  erossover  study.  The  configuration  in  Case  D2  was  selected  for  use  in  the  next  test 
series  investigating  the  shock-refleeting  surface  because  the  instantaneous  expansion  required  the 
least  downstream  distance  to  implement  and  will  be  shown  to  have  the  largest  separation 
between  decoupled  shock  and  flame  fronts  believed  to  be  important  to  prevent  trapping 
seeondary  detonations  between  the  reflecting  surface  and  the  flame  front. 


Table  1.  Diffraction  cases 


Case  # 

Diffraction  angle  (°) 

Corner  radius  (mm) 

D1 

0 

00 

D2 

90 

2.0 

D3 

15 

2.0 

D4 

90 

25.4 

A  seeond  set  of  eight  test  cases  examined  the  limits  of  detonation  reinitiation  caused  by 
oblique  shoek  refleetion.  The  parameters  of  the  set  were  primarily  reflection  angle  (P), 
downstream  distanee  from  the  diffraetion  eomer  (xo)  and  vertical  distance  from  the  diffraetion 
comer  (yo)  (Fig.  29).  Case  R5  revisited  diffraetion  comer  radius  (r)  to  mle  out  any  effect  on 
reinitiation  not  inferable  from  the  diffraction  cases.  Table  2  lists  the  eonfigurations  tested. 
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Cases  R1  and  R2  bounded  the  downstream  limit  for  placement  of  a  reflecting  surface.  Cases  R2, 
R3,  and  R4  found  a  minimum  limit  for  the  reflecting  surface  angle.  Cases  R3  and  R5  looked  for 
any  relationship  between  diffraction  comer  radius  and  reinitiation.  Cases  R6,  R7,  and  R8 
bounded  two  limits  on  reinitiation  for  a  reflecting  surface  offset  50.8  mm  vertically  from 
diffraction  comer.  The  reflection  cases  did  not  include  investigation  of  the  diffraction  of 
secondary  detonations  that  form  due  to  shock  reflection. 


r 


Entrance  Section 


xO  - 

Variable  Section 


Figure  29.  Reflection  case  parameters 


Table  2.  Reflection  cases 


Case 

r  (mm) 

P(°) 

xO  (mm) 

yO  (mm) 

R1 

2.0 

45 

162.2 

-50.8 

R2 

2.0 

45 

84.7 

-50.8 

R3 

2.0 

30 

80.1 

-50.8 

R4 

2.0 

15 

0.0 

-50.8 

R5 

25.4 

30 

80.1 

-50.8 

R6 

2.0 

-45 

0 

50.8 

R7 

2.0 

-45 

43.0 

50.8 

R8 

2.0 

-45 

169.3 

50.8 

A  third  set  of  11  cases  used  the  information  on  angle  and  position  to  investigate 
configurations  for  a  detonation  diffuser.  The  third  set  of  cases  followed  an  iterative  path  starting 
from  Case  R2  with  each  step  improving  the  chances  of  reinitiation  and  reducing  the  chance  of 
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decoupling.  Table  3  lists  the  configurations  and  their  design  parameters.  The  first  diffraction 
angle  and  comer  radius  were  constant  for  all  cases  at  90°  and  2.0  mm  respectively.  Case  Ml 
began  by  attempting  to  reduce  or  eliminate  decoupling  after  a  secondary  detonation  formed  by 
reducing  the  height  of  the  obstacle  (h)  and  increasing  the  radius  of  the  secondary  diffraction 
comer  (r2).  As  the  tests  cases  progressed,  others  design  variables  were  included  such  as  the 
number  of  reflecting  surfaces,  their  position,  and  the  number  steps  in  the  initial  diffraction.  In 
Table  3,  the  design  parameters  are  highlighted  when  they  change  to  emphasize  the  design 
choices. 


Table  3.  Multi-obstacle  test  cases 


Case 

Obstacle 
Height  (mm) 

1*2 

(mm) 

Obstacle 

Count 

Final 
Channel 
Height  (mm) 

Top  or 
Bottom 

Diffraction 

Steps 

Number  of 
Trials 

R2 

50.8 

0.3 

1 

241 

Bottom 

1 

Ml 

12.7 

6.4 

1 

241 

Bottom 

1 

4 

M2 

6.21 

0.3 

14 

241 

Bottom 

1 

22 

M3 

5.43 

3.2 

14 

241 

Bottom 

1 

7 

M4 

6.21 

0.3 

14 

102 

Bottom 

1 

8 

M5 

5.43  (top)  / 
6.21  (bottom) 

3.2/ 

0.3 

28 

102 

Both 

1 

16 

M6 

5.43 

3.2 

28 

102 

Both 

1 

7 

M7 

6.35 

0.3 

12 

102 

Top 

8 

10 

M8 

6.35 

0.3 

12 

102 

Bottom 

8 

5 

M9 

6.15 

0.3 

5 

241 

Both 

3 

9 

MIO 

6.15 

0.3 

4 

241 

Both 

3 

5 

Mil 

N/A 

0.3 

4 

241 

N/A 

1 

8 

Due  to  the  number  of  variables  considered  in  the  M-series  cases  and  the  iterative 
development  of  each,  it  is  useful  to  describe  the  cases  in  turn.  Case  Ml  (Fig.  30)  reduced  the 
height  of  the  obstacle  (hi)  in  Case  R2  from  50.8  mm  to  12.7  mm  to  verify  that  a  smaller  obstacle 
also  reinitiates  detonation,  and  increased  the  diffraction  comer  radius  (r2)  at  the  end  of  the 
reflecting  surface. 
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Figure  30.  Case  Ml 


The  decision  to  change  two  design  variables  in  parallel  would  seem  to  prevent  independent 
analysis  of  each  variable,  however;  the  diffraction  comer  radius  had  no  bearing  on  the  chances 
for  reinitiation  since  the  leading  shock  encountered  it  after  the  reflecting  surface.  Decoupling  of 
the  secondary  detonation  was  only  applicable  once  reinitiation  occurred,  and  the  only 
prerequisite  for  comparing  secondary  diffraction  comers  was  a  coupled  detonation  front.  As  a 
result,  the  secondary  diffractions  are  comparable  between  cases  where  the  obstacle  height 
differed.  After  successful  reinitiation  in  Case  Ml,  the  secondary  detonation  decoupled  both  at 
the  secondary  diffraction  comer  and  in  the  region  above  the  obstacle.  Case  M2  continued  to 
reduce  the  height  of  the  obstacle,  and  included  more  obstacles  in  an  attempt  to  increase  the 
chances  of  a  secondary  detonation  surviving  (Fig.  31). 
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There  were  14  separate  reflecting  surfaces  in  Case  M2.  This  configuration  was  expected 
to  cause  a  reinitiation  event  at  every  reflection  surface.  Then  the  secondary  detonations  would 
each  partially  decouple  before  the  next  reinitiation.  Each  detonation  then  followed  a  series  of 
decoupled  shocks  as  they  traversed  the  lead  shock  front.  Once  a  sufficient  number  of  transverse 
shocks  were  present,  their  combined  compression  would  allow  a  secondary  detonation  to 
propagate  completely  across  the  lead  shock  front  completing  the  transition  from  subcritical  to 
supercritical  channel  height. 

In  testing,  reinitiation  occurred  for  an  average  of  1 1  obstacles  before  the  lead  shock  lost 
too  much  strength.  The  compression  of  the  shocks  resulting  from  the  preceding  secondary 
detonations  was  insufficient  for  the  last  one  to  remain  coupled,  and  detonation  did  not  propagate 
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downstream.  To  mitigate  some  of  the  decay  of  the  secondary  detonations  the  tops  of  the 
obstacles  were  rounded  to  increase  the  secondary  diffraction  radii  in  case  M3. 

For  M3,  increasing  the  secondary  diffraction  comer  radii  from  0.3  mm  to  3.2  mm 
increased  the  mean  number  of  reinitiations  from  10  to  11.  While  the  rj  change  was  a  move  in  the 
right  direction,  it  was  far  from  enough  improvement  for  detonation  to  survive.  Because  the  lead 
shock  reached  the  end  of  the  obstacles  before  the  first  transverse  shock  reached  the  top  of  the 
charmel,  the  author  hypothesized  that  reducing  the  channel  height  (M4)  would  allow  the  later 
transverse  shocks  to  interact  with  the  reflections  of  the  preceding  shocks  before  the  end  of  the 
test  section.  Further,  reducing  the  final  channel  height  (M4)  decreased  the  decay  of  the 
transverse  shocks  since  the  shocks  traversed  the  lead  shock  in  less  time.  It  was  hypothesized  that 
the  shock  collisions  and  reduced  traverse  time  would  cause  local  explosions  along  the  top  of  the 
charmel,  and  that  the  secondary  detonations  from  the  top  and  bottom  of  the  channel  would  merge 
into  a  single,  fully  coupled  detonation  front.  Testing  in  Case  M4  did  not  confirm  the  hypothesis 
because  no  local  explosions  were  observed  at  the  top  of  the  charmel  likely  because  the  shock 
decay  was  still  sufficient  to  reduce  the  probability  of  reinitiation  to  zero.  Another  iteration  was 
necessary  and  Case  M5  added  obstacles  to  the  top  of  the  channel. 

For  M5,  because  obstacles  reliably  reinitiated  detonation  from  a  decoupled  shock  at  the 
bottom  of  the  charmel  when  the  first  obstacle  was  85  mm  downstream  of  the  diffraction  comer,  it 
was  reasonable  to  assume  the  same  was  tme  of  the  top  of  the  channel  when  the  first  obstacle  was 
closer  to  the  diffraction  comer.  A  new  set  of  obstacles  identical  to  that  of  case  M3  were  placed 
at  the  top  of  the  channel  without  changing  the  final  charmel  height  from  M4  (Fig.  32).  This 
configuration  was  the  first  multi-obstacle  configuration  to  exhibit  shock  initiated  combustion 
along  the  top  of  the  channel. 
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yO  =  50.8  mm 

^ ^ ^ 

Obstacle  Set  from  Case  M3 

xO  -  0.00  mm  i 
yO  =  -50.8  mm 

Obstacle  Set  from  Case  M2 

^ 

Entrance  Section 

_ Variable  Section _ 

Figure  32.  Case  M5 


In  one  of  the  eight  repeat  trials  of  the  M5  configuration,  the  1 1*  obstacle  reinitiated  a  detonation 
that  failed  to  reach  the  end  of  the  test  section.  The  results  in  Case  M5  were  promising,  and  it 
seemed  reasonable  to  round  the  diffraction  comers  of  the  obstacles  at  the  bottom  of  the  channel 
again  for  Case  M6. 

As  in  Case  M3,  increasing  the  diffraction  radii  did  not  cause  enough  improvement  for 
successful  frilly  coupled  detonation  at  the  end  of  the  test  section.  One  of  the  16  repeat  trials  for 
M6  resulted  in  a  secondary  detonation  at  the  top  of  the  channel.  Unlike  in  case  M5,  the 
detonation  remained  partially  coupled  through  the  end  of  the  test  section.  At  this  point,  it 
seemed  unlikely  that  obstacles  set  so  far  vertically  from  the  initial  diffraction  comer  would  be 
sufficient  to  reinitiate  detonation  by  the  end  of  the  test  section.  Case  M7  began  to  modify  the 
diffraction  comer  by  dividing  the  change  in  channel  height  into  eight  discrete  segments  (Fig.  33). 
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the  8.19  mm  cell  size  of  stoichiometric  hydrogen/air  at  1  atm  initial  pressure  (Ciccarelli,  1994). 
Setting  the  step  height  equal  to  the  cell  width  allowed  an  average  of  one  triple  point  to  interact  at 
each  step.  The  reflected  triple  points  were  expected  to  maintain  detonation  unlike  in  previous 
cases.  A  new  set  of  12  obstacles  followed  the  diffraction  steps  along  the  top  of  the  channel.  The 
obstacles  had  a  shorter  pitch  of  12.4  mm  to  increase  the  rate  of  diffraction  and  reinitiation.  In 
testing,  the  primary  detonation  decoupled  after  the  first  diffraction  step  and  shock  reflections 
from  the  obstacles  failed  to  reinitiate  detonation.  Since  the  obstacles  were  unable  to  reinitiate 
detonation  at  the  top  of  the  channel,  they  were  moved  to  the  bottom  of  the  channel  in  Case  M8 
where  it  was  known  that  reinitiation  would  occur. 
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The  obstacles  reinitiate  detonation  in  Case  M8;  however,  the  secondary  detonations  again 
decoupled  before  reaching  the  top  of  the  channel.  The  decoupling  along  the  diffraction  steps 
contradicted  the  hypothesis  that  the  short,  cell-sized  steps  would  partially  mitigate  decoupling  of 
the  primary  detonation  front.  Instead  of  mitigating  the  decoupling,  it  was  reasoned  that  a  new 
diffraction  wall  (Fig.  34)  could  be  redesigned  to  reflect  the  decoupled  lead  shock  twice  for  each 
diffraction  comer  as  shown  by  the  outline  in  Figure  30.  Case  M9  implemented  a  double 
reflection  both  on  the  diffraction  wall  and  on  the  bottom  of  the  channel  (Fig.  35). 


Diffraction  Wall 


Entrance  Section 


Variable  Section 


Figure  34.  Diffraction  wall  location 
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Early  in  the  testing  of  Case  M9,  it  became  apparent  that  the  first  obstacle  on  the  bottom 
of  the  channel  was  restricting  the  flow  of  reactants  in  to  the  test  section.  As  a  result,  the  initial 
detonation  decoupled  before  reaching  the  test  section.  To  reduce  the  restriction,  the  first  obstacle 
on  the  bottom  of  the  channel  was  removed  for  Case  MIO.  Removing  the  obstacle  reduced  the 
restriction  enough  that  the  initial  detonation  no  longer  decoupled  before  the  test  section. 

In  testing  MIO  detonation  reinitiation  occurred  at  the  first  bottom  obstacle  in  each  of  five 
repeat  trials,  but  nowhere  else.  There  was  some  infrequent  shock  ignition  along  the  upper 
obstacles.  At  this  point  lining  the  upper  and  lower  walls  of  the  channel  with  obstacles  had  failed 
to  do  more  than  reinitiate  one  secondary  detonation  at  a  time  none  of  which  were  sufficiently 
strong  to  complete  a  transition  to  detonation  in  the  final  channel.  The  next  logical  step  was  to 
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reinitiate  several  detonations  at  roughly  the  same  time  and  then  let  them  merge  into  a  single 
detonation  front. 

In  Case  Ml  1  (Fig.  36),  four  obstaeles  were  attaehed  to  the  windows  of  the  test  seetion  so 
that  the  channel  split  into  five  channels.  The  obstacles  were  arranged  using  the  shock  Mach 
number  map  from  case  D2  so  that  reinitiation  would  occur  on  one  wall  of  each  of  the  five 
channels.  Then  the  secondary  detonations  would  travel  down  the  five  channels  expanding 
gradually  before  emerging  at  the  same  time  to  merge  into  a  single  detonation  front. 


This  arrangement  resulted  in  the  first  local  explosions  to  occur  as  the  result  of  collisions  of  two 
diffracting  shockwaves  (see  Fig.  78).  The  behavior  appears  identical  to  that  of  detonations 
diffracting  after  emerging  from  critical  channels  (Soloukhin  and  Ragland,  1965).  The  Mil 
configuration  bridged  the  gap  between  the  subcritical  and  critical  cases  of  detonation  diffraction. 
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In  Chapter  III,  the  experimental  hardware  was  discussed  at  some  length.  The  focus  was 
on  the  methods  for  obtaining  a  repeatable  experiment  and  providing  complete  descriptions  of  the 
test  cases  and  the  reasoning  behind  the  design  choices  in  each  case.  In  Chapter  IV  the  focus 
shifts  to  the  collection  and  analysis  of  data  collected.  The  image  analysis,  data  reduction,  and 
uncertainty  calculations  give  quantitative  meaning  to  the  results  presented  in  Chapter  V. 
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IV.  Analysis  Methodology 


Overview 

Schlieren  images  of  shocks  and  flames  were  used  to  construct  a  clear  picture  of  the 
diffraction  and  reinitiation  processes.  Tedious  hand  selection  of  points  along  the  shock  and 
flame  fronts  gave  a  basis  for  analysis  of  the  motion  of  those  structures.  The  distance  between  the 
flame  and  shock  were  derived  from  a  forward,  closest-point  measurement  (see  Fig.  9)  from  the 
flame  to  the  shock.  The  velocity  of  the  shock  front  was  derived  from  an  interpolation  of  the 
shock  position  in  two  adjacent  frames.  These  two  measurements  were  sufficient  to  construct 
accurate  interpolation  functions  that  enabled  the  prediction  of  reinitiation. 

Shock  and  Flame  Position 

The  first  step  in  the  analysis  was  hand  selection  of  the  shock  and  flame  position.  Shocks 
were  identifiable  by  a  thin,  smooth  wave  front  at  a  position  ahead  of  any  flames  (Fig.  37). 


Shock 


Flame 


Figure  37.  Shock  and  flame  fronts 


52 


Flames  were  much  thicker  than  shocks  and  usually  included  protrusions  resulting  from  local 
variations  in  flame  speed.  Efforts  to  automate  the  identification  of  shocks  and  flames  in  the 
schlieren  images  were  thwarted  by  the  presence  of  reflected  shocks,  detonation  fronts  and 
multiple  flame  fronts,  and  the  author  resorted  to  hand  processing  of  the  schlieren  images.  An 
author-built  software  tool,  included  in  Appendix  C  greatly  aided  the  process.  The  tool  loaded  the 
desired  frame  from  a  video  and  modified  it  for  display.  The  modification  process  subtracts  a 
background  image  of  window  defects,  adjusts  the  contrast,  and  shades  the  visible  walls  and 
obstacles  blue.  The  tool  then  displays  the  modified  frame  (Fig.  38)  and  prompts  the  user  to 
select  the  first  point  on  either  a  shock  or  flame.  Before  making  a  selection,  the  user  has  the 
option  to  change  the  magnification  and  position  of  the  window  in  relation  to  the  image  for  ease 
of  visibility. 


Figure  38.  Initial  state  of  software  tool  interface  showing  diffraction  corner  and  enhanced  image 
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Once  satisfied  with  the  window  size  and  position,  the  user  selects  the  first  pixel  on  one 
end  of  the  shock  or  flame.  The  software  colors  that  pixel  red  (Fig.  39).  Then,  it  records  the 
coordinates,  and  changes  the  input  mode  from  mouse  to  keypad. 


Figure  39.  First  point  on  shock  selected,  ready  to  select  adjacent  pixel. 

In  keypad  mode,  pressing  a  key  adds  the  coordinates  of  the  adjacent  pixel  in  that  direction  to  the 
list  of  coordinates  describing  the  shock  or  flame.  Pressing  "4"  adds  the  pixel  directly  to  the  left 
and  pressing  "6"  adds  the  pixel  to  the  right. 

Pixel  selection  continues  as  the  user  selects  the  next  pixel  along  the  structure.  With  each 
selection,  the  software  marks  the  new  pixel,  updates  the  list  of  coordinates,  and,  when  necessary, 
re-centers  the  view.  Mistakes  due  to  typos  (Fig.  40)  were  corrected  via  an  undo  feature  that 
rewinds  the  selection  one  pixel  at  a  time  by  pressing  the  5  key. 


54 


Figure  40.  Mistakenly  pressing  9  instead  of  6  results  in  this  state. 

The  result  is  a  path  of  red  pixels  on  screen  and  an  ordered  list  of  the  coordinates  of  the 

pixels  (Fig  41).  When  the  user  completes  the  selection,  the  software  returns  the  list  of 
coordinates  and  exits.  Flames  were  handled  in  the  same  manner,  by  following  the  leading  edge 
of  the  combustion  front  including  any  protrusions  and  recesses.  Runtime  per  frame  is  around  2 
minutes,  and  the  2002  frames  processed  in  this  manner  yielded  about  509,000  coordinate  pairs. 
The  tool  took  about  10  hours  to  code,  and  the  total  time  to  hand  process  the  shock  and  flame 
positions  was  approximately  180  man-hours. 


Figure  41.  Completed  selection  of  the  shock  front  and  the  first  flfteen  coordinate  pairs 
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Separation  Distance 

A  forward,  closest-point  method  calculated  the  separation  distance  for  each  point  on  the 
flame  front.  The  algorithm  first  applied  a  calibration  (Eq.  4)  to  the  shock  and  flame  pixel 
coordinates  to  convert  from  pixels  to  millimeters. 

_  ^initial,  mm 

•^mm  y.  -^pix  V  V 

'Hnitial,  pix 

The  algorithm  then  calculated  the  distance  from  each  point  on  the  flame  to  all  of  the  points  on 
the  shock  using  Eq.  5. 


distance  =  +  Ay^  (5) 

The  algorithm  then  assumes  that  the  minimum  distance  best  represents  the  separation  distance  at 
that  location  on  the  flame  (Fig.  42). 


(mm) 

Figure  42.  Separation  distance  vectors 

The  bias  uncertainty  in  Ax  and  Ay  is  half  the  pixel  width  (0.324  mm)  and  the  bias  uncertainty  in 
distance  is  given  in  Eq.  6. 
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^distance  =  {{Bax  '  Ax/^Ax^  +  Ay^)  +  {B^y  ■  Ay/^Ax^  +  Ay^)  )  ^  (6) 

for  a  sample  where  Ax  =  3  pix  and  Ay  =  4  pix,  the  separation  distance  is  5.00  ±  0.324  mm. 

Shock  Mach  Number 

The  post-shock  flow  field  (pressure,  temperature,  density,  and  velocity)  depends  solely 
on  the  Mach  number  at  which  the  wave  travels.  Mach  number  was  calculated  from  the 
temperature  of  the  reactants,  the  equivalence  ratio,  stoichiometry,  and  the  frame -to-frame  shock 
displacement  as  shown  below.  The  post-shock  conditions  were  calculated  from  the  shock  jump 
equations  (Anderson,  1982). 

The  sole  mixture  used  in  experiments  was  stoichiometric  hydrogen  and  air.  At 
stoichiometric  conditions,  the  equivalence  ratio  ((|))  is  unity.  For  each  mole  of  hydrogen 
combusted,  the  equivalence  ratio  dictates  the  number  of  moles  of  oxygen  and  nitrogen  present 
(Eq.  7)  (Turns,  2000). 

No2  =  (x  +  yl4)lcl>  (7) 

where  x  is  the  number  of  carbon  atoms  (0)  and  y  is  the  number  of  hydrogen  atoms  (2)  per 
molecule  of  fuel.  The  number  of  moles  of  nitrogen  is  determined  by  the  natural  ratio  of  nitrogen 
to  oxygen  in  air  (Eq.  8). 

N^2  =  3.76  *  No2  (8) 

The  mole  fractions  of  the  various  reactants  are  equal  to  the  number  of  moles  of  that  species 
divided  by  the  total  number  of  moles  (Eq.  9). 

Xi  =  Ni/T,iNi  (9) 
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The  bias  uncertainties  in  the  Ni  and  Xi  were  calculated  using  Eq.  10  (Coleman  and  Steele,  1989). 


Bf]  — 


Table  4  gives  values  for  Ni  and  Xi  and  their  respective  uncertainties. 


(10) 


Table  4.  Stoichiometry  variable  sample  values 


Variable 

value 

uncertainty 

Nh2 

1.0 

0  (exact) 

1.000 

0.001 

No2 

0.5 

0.001 

Nn2 

1.88 

0.00188 

ENi 

3.38 

0.00195 

XH2 

0.2959 

1.604E-4 

X02 

0.1479 

1.539E-4 

XN2 

0.5562 

5.788E-4 

The  specific  heat  at  constant  pressure  of  the  mixture  is  a  function  of  the  mole  fractions 
and  the  specific  heats  of  the  component  species  (Eq.  1 1). 

^p,mix  ~  (H) 

The  specific  gas  constant  for  the  mixture  is  a  function  of  the  universal  gas  constant,  the  mole 
fractions,  and  the  molecular  weights  of  the  components  (Eq.  12). 

R  =  Ru/I.iMWi-Xi  (12) 

The  ratio  of  specific  heats  is  a  function  of  Cp  and  R  (Eq.  13) 

7  =  Cp/(Cp-i?)  (13) 

The  speed  of  sound  in  the  test  section  (a)  is  a  function  of  the  ratio  of  specific  heats  (y), 
temperature  (T),  and  the  mixture  specific  gas  constant  (R)  (Eq.  14).  The  speed  of  sound  was 
constant  throughout  the  test  cases.  Table  5  gives  the  values  and  uncertainties  for  Cp,  R,  y,  and  a 
in  addition  to  the  MWi  and  Cp,i  values  used  in  Eqs.  1 1  and  12. 

a  =  ^y-R-T  (14) 
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Table  5.  Sample  Mach  number  and  constituent  values 


Variable 

Value 

Bias  Uncertainty 

Units 

Cp,H2 

28.877 

0.001 

kJ/kmol-K 

Cp,02 

29.331 

0.001 

kJ/kmol-K 

Cp,N2 

29.075 

0.001 

kJ/kmol-K 

MWh2 

2.016 

0.001 

g/mol 

MWo2 

31.999 

0.001 

g/mol 

MWn2 

28.013 

0.001 

g/mol 

Cp 

1.3811 

0.0156 

kJ/kg-K 

R 

0.3976 

0.00357 

kJ/kg-K 

y 

1.4043 

0.0384 

- 

a 

409.3 

0.144 

m/s 

The  Mach  number  of  the  shock  was  found  by  dividing  the  distance  traveled  between 
frames  by  the  time  interval  and  the  speed  of  sound  (Eq.  15). 

M  =  +  Ay^/ (At  •  a)  (15) 

The  distance  traveled  by  the  shock  was  found  using  the  central  finite  difference  method 
described  in  Figure  43. 


1.  Identify  shock  location  in  frame  . 

(redj  6  squares) 

2.  Identify  shock  position  in  2^^  frame  . 
(blue,  9  squares) 

3.  Interpolate  to  acquire  an  equal  number 
of  points  for  each  shock,  (dots) 

4.  Average  the  interpolated  positions  to  get 
average  position  during  interval. 

(black  dots) 

5.  Measure  displacement  (black  arrows)  and 
divide  by  At  to  get  shock  speed. 


Figure  43.  Algorithm  for  determining  shock  speed. 
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The  central  difference  method  reports  speeds  half  way  between  the  two  shocks  and  requires  an 
equal  number  of  starting  and  ending  points.  Unlike  flames,  shocks  have  a  relatively  slowly 
varying  curvature,  and  it  was  appropriate  to  interpolate  points  along  the  shocks  so  that  there  were 
an  equal  number  of  points  on  each  shock  front.  Averaging  the  locations  of  the  two  shocks 
yielded  the  midpoints  where  the  Mach  number  was  reported.  The  bias  uncertainty  in  the 
horizontal  and  vertical  distances  was  the  width  of  a  pixel,  and  the  manufacturer-reported 
uncertainty  in  frame  interval  (At)  was  20  ns.  Table  6  reports  sample  values  and  bias 
uncertainties  of  Ax,  Ay,  At  and  M. 


Table  6.  Sample  Mach  number  and  constituent  variables 


Variable 

Value 

Bias  Uncertainty 

Units 

a 

409.3 

0.144 

m/s 

Ax 

7.958*10'^ 

0.647*10'^ 

m 

Ay 

6.580*10’^ 

0.647*10’^ 

m 

At 

4.75*10'^ 

20*10'^ 

s 

M 

5.312 

0.334 

Interpolation  Functions 

Measurements  of  flame  separation  and  shock  Mach  number  could  not  be  collected  for  a 
regularly  spaced  grid  of  positions  due  to  the  frame  rate  limit  for  the  Phantom  v71 1  camera.  The 
time  interval  between  frames  resulted  in  discrete  spatial  measurements.  As  a  result,  there  were 
spatial  gaps  in  the  data  (Fig.  44),  and  very  few  locations  had  the  repeated  measurements  needed 
to  calculate  precision  error.  To  reconcile  the  statistical  variation  and  fdl  in  the  voids,  a  fitting 
function  was  devised.  The  fitting  function  had  to  be  continuous  in  x  and  y  and  be  a  linear 
combination  of  terms  such  that  a  linear  least  squares  fit  could  be  applied  to  find  the  coefficients 
for  each  term.  A  two-dimensional  power  series  (Eq.  16)  met  the  requirements,  and  linear  least 


60 


squares  fitting  was  used  to  obtain  eoeffieients  (an,m)  from  the  measured  separation  distanee  and 
Maeh  number  data  in  eaeh  test  ease. 


mm 

Figure  44.  Consolidated  shock  position  measurements  from  all  runs  of  Case  D2.  Note  spatial  gaps  in 
measurements  due  to  40  kfps  frame  rate  used  to  increase  image  size. 

f{x.  y)  =  i:n=0  lm=0  an,mX^y^  (16) 

In  practice,  the  series  must  be  truncated  to  a  maximum  power  (p).  The  choice  of  p  is  a  tradeoff 
between  minimizing  the  error  of  the  function  at  the  measured  points  and  limiting  the  effect  of 
outliers  on  the  function.  In  each  test  case,  the  series  was  truncated  to  the  power  that  produced 
the  least  mean  absolute  error. 

The  coefficients  were  calculated  using  the  linear  least  squares  regression.  The  regression 
seeks  to  minimize  the  difference  between  the  fitting  function  and  the  measured  data.  It  takes  the 
form  of  Eq.  17 
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y=[X]a  +  8  (17) 

where  the  elements  of  y  are  the  measured  separation  distances  or  Mach  numbers,  a  is  a  vector  of 
the  coefficients,  and  £  is  the  error.  For  a  data  set  with  n  measurements  and  a  fitting  function  of 
power  p,  [X]  is: 

1  •••  yf' 

[X]  =  ^  (Ig) 

-1  •••  Vn- 

Equation  19  gives  the  coefficients  that  minimize  the  root  sum  square  of  the  elements  of  £. 

a=[X]~'^*y  (19) 

The  full  code  of  the  fitting  function  is  included  in  Appendix  C. 

The  fitting  functions  allowed  the  calculation  of  precision  errors  based  on  the  entire  data 
set  for  each  test  case.  When  calculating  the  precision  error,  the  fitted  data,  Xfn,  was  substituted 
for  the  mean,  X,  in  Coleman  and  Steele  to  give  Equation  20  (Coleman  and  Steele,  1989).  The 
value  of  t  in  Eq.  20  was  1.96  for  95%  a  confidence  interval  because  all  of  the  data  sets  had 
sufficient  samples  to  use  the  normal  distribution  in  place  of  the  student's  t-distribution. 

Px  =  t*  \  /VH  (20) 

The  maximum  precision  error  in  separation  distance  was  6.64x10'  mm  and  the  maximum 
precision  error  in  Mach  number  was  0.0124.  Both  of  the  precision  errors  were  an  order  of 
magnitude  smaller  than  the  bias  errors;  therefore,  the  bias  dominates  the  total  error.  Combining 
the  bias  and  precision  errors  gives  maximum  uncertainties  of  0.650  mm  (4.34%)  for  separation 
distance  and  0.334  (6.29%)  for  Mach  number. 
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V.  Results 


Overview 

A  proof  of  concept  experiment  was  designed  by  the  author  to  prove  the  existence  of 
secondary  detonations  caused  by  oblique  reflection  of  a  diffraction  detonation  (Fig.  45).  An 
increase  in  the  chemiluminescence  along  the  reflecting  surface  indicated  detonation  in  the  test 
section,  and  wave  speed  measurements  downstream  of  the  reflecting  surface  also  indicated  a 
detonation,  but  the  integration  time  of  the  images  in  this  initial  look  was  too  long  to  quantify  the 
conditions  leading  to  reinitiation. 


After  the  proof  of  concept  experiment,  the  author  began  an  organized  exploration  of 
diffraction  parameters  with  the  D-series  test  cases.  Study  of  the  diffraction  angle,  0,  and  the 
comer  radius,  r  (Fig.  46)  revealed  that  radius  had  little  effect  on  the  separation  distance  or  shock 
Mach  number  while  there  was  no  need  for  reinitiation  hardware  when  the  diffraction  angle  was 
less  than  15°  because  the  detonation  did  not  fully  decouple. 

In  the  R-series  cases,  it  was  observed  that  reinitiation  occurred  within  a  bounded  set  of 
the  reflection  ramp  parameters.  The  maximum  downstream  distance  (xo)  for  reinitiation  from  a 
45°  surface  was  between  1.67  and  3.19  times  the  initial  channel  height,  hintiai,  when  yo  =  0. 
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Figure  46.  Parameters  of  the  D-series  and  R-series  cases 


The  minimum  reflecting  angle  (P)  for  reinitiation  at  the  plane  of  the  diffraction  comer  was  15°. 
At  a  vertical  offset  (yo)  of  hintiai,  a  45°  reflecting  surface  caused  reinitiation  only  at  downstream 
distances  between  0.846  hintiai  and  3.33  hintiai.  Variation  in  the  diffraction  comer  radius  had  no 
significant  effect  on  diffraction  or  reinitiation  contrary  to  the  qualitative  evidence  in  Nielsen  et 
al.  (2011). 

Diffraction  of  the  secondary  detonations  prevented  any  of  eleven  M-series  cases  from 
successfully  propagating  detonation  into  the  final  channel  height  (Fig.  47).  All  of  the  M-series 
cases  were  analyzed  qualitatively  foregoing  a  time  consuming  quantitative  analysis  until  a 
successful  geometry  was  observed.  The  Cases  Ml  and  M2  showed  that  reinitiation  occurred  for 
reflecting  obstacle  heights  as  small  as  0.77X,  but  rounding  the  diffraction  comer  at  the  end  (Ml) 
did  not  prevent  decoupling.  Cases  M2  through  MIO  added  more  reflecting  surfaces  along  the 
walls  of  the  channel  and  those  suffered  from  a  combination  of  repeated  diffraction  of  the 
reinitiated  detonation  waves.  In  some,  more  than  one  reinitiation  occurred  but  the  detonation 
never  survived  into  the  supercritical  channel.  Interestingly,  case  Mil  displayed  a  second  round 
of  reinitiation  after  diffraction  of  the  secondary  detonations  (Fig.  48).  In  this  case,  reinitiation 
was  caused  by  the  collision  of  two  decoupled  shocks  from  two  secondary  detonations  in  adjacent 
radial  channels. 
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Figure  47.  Diffraction  of  the  secondary  detonation  causes  it  to  decouple 
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Diverging  and  Converging  Channel  Tests 

Stevens  et  al.  (2011)  eontains  the  first  experimental  data  eollected  specifically  for  the 
development  of  a  detonation  diffuser.  Two  geometric  test  cases  were  considered,  a  14° 
diverging  ramp  and  a  step  expansion  followed  by  a  14°  converging  ramp  (Fig.  49).  The  two 
geometries  utilize  different  approaches  to  a  diffuser.  The  diverging  case  seeks  to  limit  the 
diffraction  angle  such  that  the  initial  detonation  never  decouples,  and  the  converging  case  allows 
decoupling  to  occur  so  that  detonation  will  reinitiate  due  to  the  shock  reflection  from  the 
converging  wall. 


Since  the  cell  size  was  8.19  mm,  the  initial  channel  height  of  38.1  mm  was  subcritical 
(h/A,  =  4.7)  in  both  cases.  The  diffraction  angle  in  the  diverging  case  (14°)  was  0.5°  less  than  the 
predicted  limit  of  14.5°  (Nettleton  1987).  The  final  channel  height  of  50.8  mm  was  barely 
supercritical  (h/A  =  6.20)  according  to  Lee’s  (1995)  definition  of  the  critical  height,  but  it  was  as 
large  as  the  test  section  allowed  at  the  time.  The  length  of  the  transition  to  a  supercritical 
channel  would  be  176  mm. 
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Schlieren  imaging  of  the  diverging  case  indicated  decoupling  despite  the  small  ramp 
angle  (Fig.  50).  To  avoid  decoupling,  the  diffraction  angle  has  to  be  smaller,  but  other  angles 
were  not  tested.  For  a  benchmark  to  compare  to  the  converging  case,  it  was  assumed  that  a  13.5° 
diffraction  angle  was  sufficient  to  avoid  decoupling  in  the  diverging  case  since  14°  was  not 
sufficient.  The  target  transition  length  was  182  mm  corresponding  to  the  13.5°  diffraction  angle. 
It  was  hypothesized  that  a  diffract,  decouple,  and  reinitiate  approach  would  complete  the 
transition  in  less  distance  due  to  the  shorter  distance  needed  for  a  90°  diffraction.  This  reasoning 
lead  to  the  exclusive  use  of  diffraction  and  reinitiation  rather  than  diffraction  alone  in  the 
remainder  of  this  work. 
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Frame  2: 
13.3  |r 


Frame  3: 
26.7  |rs 


Frame  4: 
40.0  |as 


Diffracting  front 
Head  of  expansion 

Unaffected  Detonation 


Possible  local  explosion  at 
window 


Visible  decoupling 
Expansion  head 

ontinued  Decoupling 


Figure  50.  Decoupling  on  the  diverging  ramp  (Stevens  et  al.,  2011) 


In  the  converging  case,  the  initial  channel  height  was  also  38.1  mm.  The  step  widened 
the  channel  by  12.7  mm  to  its  maximum  (50.8  mm).  Then  the  converging  wall  reduced  the 
channel  height  back  down  to  38.1  mm  with  a  reflecting  angle  of  14°.  It  was  unknown  if,  or 
where,  detonation  would  reinitiate,  but  any  observed  reinitiations  had  to  occur  at  a  location 
where  the  channel  height  was  greater  than  the  initial  height. 
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Schlieren  images  indicated  reinitiation  at  the  midpoint  of  the  converging  wall  (Fig.  51). 
The  first  frame  of  the  video  shows  decoupling  of  the  primary  detonation  after  the  12.7  mm  step. 
The  second  frame  shows  two  bright  regions  where  reinitiation  occurred.  The  first  is  in  the  comer 
where  the  converging  ramp  meets  the  top  of  the  test  section.  The  second  begins  at  the  midpoint 
of  the  converging  wall  and  continues  off  to  the  right  side  of  the  frame.  The  channel  height  at  the 
point  of  reinitiation  was  44.5  mm.  The  successful  reinitiation  lead  to  some  speculation  on  the 
result  of  adding  more  obstacles  to  the  geometry.  Had  the  converging  wall  ended  at  the  point  of 
reinitiation  and  been  followed  by  another  12.7  mm  step  and  converging  wall  the  pattern  of 
diffract,  decouple,  and  reinitiate  could  have  repeated  as  many  times  as  needed  to  reach  the 
critical  channel  height.  The  total  length  of  such  a  transition  is  178  mm,  a  2.3%  reduction  in 
length.  Due  to  the  length  reduction,  diverging  geometries  were  abandoned  for  the  rest  of  the  test 
cases. 


Frame  1:  13.3  ps 


Frame  2:  26.7  ps 


Detonation  continues 
downstream 


Figure  51.  Local  explosion  on  the  converging  ramp  (Stevens  et  al.,  2001) 
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The  diverging  and  eonverging  eases  were  useful  for  deciding  to  further  study  diverging 
or  converging  geometries,  but  the  frame  rate  was  too  low  and  the  exposure  too  long  to  accurately 
measure  the  position  of  shocks,  flames,  and  detonations.  Beginning  with  the  diffraction  series, 
the  exposure  of  the  schlieren  images  was  short  enough  and  the  spatial  resolution  (mm/pixel) 
large  enough  that  the  images  could  be  analyzed  quantitatively  with  uncertainties  sufficiently 
small  for  statistical  significance.  Due  to  the  quantity  of  data,  the  raw  images  of  the  D,  R,  and  M 
series  runs  are  included  in  Appendix  B. 

Diffraction  angle  and  diffraction  corner  radius,  cases  D1-D4 

The  D-series  of  high  precision,  quantitative  test  cases  studied  two  geometric  parameters 
of  the  diffraction  comer,  the  diffraction  angle  and  comer  radius.  Diffraction  angle  was  included 
because  it  can  unilaterally  determine  if  decoupling  occurs.  Comer  radius  was  included  because 
it  influenced  the  separation  distance  in  crossover  tubes.  These  two  parameters  were  studied  first 
because  they  do  not  require  a  reflecting  surface  in  the  test  section. 

Case  D1 

Case  D1  was  the  control  for  the  quantitative  test  cases,  a  straight  channel  (Fig.  52a). 

A  straight  channel  should  not  exhibit  diffraction  or  decoupling  of  the  detonation  front,  and  the 
wave  speed  through  the  section  is  the  Chapman-Jouguet  speed  for  the  mixture.  A  declining 
wave  speed  or  any  separation  indicates  a  poor  transition  from  the  upstream  detonation  tube  to  the 
test  section.  Figure  52b  shows  the  separation  distance  profile  for  a  combination  of  four  mns. 
One  mn  had  a  small  region  where  the  upper  part  of  the  wave  was  initially  decoupled,  visible  in 
the  upper  left  of  the  figure.  In  the  black  regions  of  the  figure,  excess  sealant  or  the  walls  of  the 
channel  blocked  the  schlieren  light  path.  The  gray  regions  were  visible,  but  fell  outside  of  the 
region  bounded  by  measured  data  points.  Outside  the  boundary,  the  fitting  function  discussed  in 
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Chapter  IV  would  rely  on  extrapolation  instead  of  interpolation  removing  the  bounds  on 
precision  error.  The  wave  quickly  recoupled,  but  it  increased  the  fitted  separation  and  lowered 
the  fitted  Mach  number  in  the  upper  left  comer  of  the  fitted  data  (Figs.  52b  and  52c).  The  mean 


Figure  52.  Case  D1  (0  =  0°,  r  =  oo)  schematic  and  data  fits  from  7  runs. 
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Case  D2 


Case  D2  had  a  diffraction  angle  of  90.0°  and  a  comer  radius  of  2.00  mm  (Fig.  53a).  This 
is  a  frequent  geometry  in  diffraction  studies  (Mitrofanov  1965,  Moen  et  al.  1982,  Shepherd  and 
Akbar  1999,  Pintgen  2004,  Arienti  and  Shepherd  2005,  Nielsen  2011,  Camardo  2012),  but  this  is 
the  first  study  to  record  the  separation  distance  in  addition  to  Mach  number. 

The  detonation  quickly  began  to  decouple  after  the  diffraction  comer  (Fig.  53b)  and  the 
entire  front  decoupled  within  30  mm  along  the  vertical  wall.  Coupled  combustion  endured 
longer  on  the  horizontal,  straight  wall  owing  to  the  sonic  propagation  of  the  expansion  across  the 
detonation  front,  but  full  decoupling  still  occurred  by  the  time  the  wave  traveled  180  mm 
downstream. 

Mach  number  degraded  quickly  along  the  diffraction  wall  as  the  separation  increased  (Fig 
53c).  This  configuration  favors  placing  a  reflecting  surface  directly  in  the  path  of  the  emerging 
detonation  front  for  reinitiation.  The  Mach  number  is  higher  along  the  bottom  wall  and  a  surface 
attached  to  the  bottom  of  the  channel  cannot  trap  reactants.  Along  the  diffraction  wall,  the  shock 
weakens  more  quickly,  and  an  offset  surface  above  the  diffraction  comer  is  prone  to  trapping 
reactants  due  to  the  nearly  vertical  propagation  of  the  shock.  Trapped  secondary  detonations 
appear  in  two  of  the  R-series  test  cases  both  of  which  had  their  reflecting  surfaces  attached  to  the 
top  wall  of  the  channel. 
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Figure  53.  Case  D2  (0  =  90°,  r  =  2  mm)  schematic  and  data  fits  from  10  runs 
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In  case  D3,  the  diffraction  angle  was  15°,  and  the  comer  radius  was  2  mm  (Fig.  54a).  As 
shown  in  Figure  54b,  the  detonation  never  fully  decouples.  Figure  54b  shows  that  some 
decoupling  occurs  because  the  separation  distance  is  greater  than  zero  in  much  of  Figure  54b,  but 
the  shock  and  flame  never  completely  separate.  The  mean  separation  remains  low,  but  is  higher 
than  in  the  straight  channel.  Since  detonation  never  fully  decouples,  the  Mach  number  remains 
high  on  average  (Fig.  54c)  with  a  slow  region  just  downstream  of  the  diffraction  comer  and  a 
gradual  decrease  as  the  channel  height  increases.  The  gradual  decrease  in  Mach  number 
suggests  that,  given  enough  length,  full  decoupling  will  occur.  Obviously,  this  arrangement  is 
ideal  for  the  survival  of  the  detonation;  however,  the  increase  in  area  is  very  gradual  and  the 
mixture  must  remain  constant  through  the  transition.  Shock  reinitiation  has  neither  of  these 
constraints. 
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Figure  54.  Case  D3  (0  =  15°,  r  =  2  mm)  schematic  and  data  fits  from  4  runs 
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Case  D4 


Case  D4  was  similar  to  case  D2  except  that  the  diffraction  corner  radius  was  25.4  mm 
instead  of  2  mm  (Fig.  55a).  The  increased  comer  radius  resulted  in  some  reduction  in  the 
separation  until  the  end  of  the  curved  section  of  wall  (Fig.  55b).  Overall,  there  was  no 
significant  reduction  in  separation.  The  Mach  number  of  the  lead  shock  was  slightly  higher  than 
that  of  D2  near  the  diffraction  comer,  but  the  increase  was  less  than  the  uncertainty  (Fig.  55c).  It 
is  possible  that  further  increase  of  the  comer  radius  could  lead  to  significant  improvement  in  the 
separation  or  the  Mach  number,  but  a  large  corner  radius  also  requires  more  fuel  to  fill  the 
volume.  For  this  reason,  comer  radius  should  not  drive  the  design  of  a  detonation  diffuser. 
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Figure  55.  Case  D4  (0  =  90°,  r  =  25.4  mm)  data  fitted  to  10  runs. 


77 


Composite  results  of  the  D-series 

In  the  diffraction  test  cases,  the  shock-flame  separation  distance  and  the  shock  Mach 
number  varied  separately.  In  cases  Dl,  D2,  and  D3  diffraction  angle  was  the  sole  independent 
variable.  Figure  56  shows  the  effect  of  diffraction  angle  on  separation  distance  (Fig.  56a)  and 
Mach  number  (Fig.  56b)  at  a  single  spatial  point  along  the  bottom  wall  of  the  test  section  200 
mm  downstream  of  the  diffraction  comer.  The  separation  distance  increases  linearly  with 
diffraction  angle  at  a  rate  of  0.1075  mm/deg.  The  correlation  coefficient  (R  )  of  the  separation 
data  is  unity  to  the  second  decimal  place,  which  indicates  that  the  trend  is  likely  linear.  The 
Mach  number  declines  with  increasing  diffraction  angle,  but  with  an  R  of  0.83,  the  trend  is  non¬ 
linear,  and  linear  trend- line  passes  outside  the  error  bars  of  two  of  the  three  points  in  Fig.  56b. 
The  decline  in  Mach  number  is  about  the  same  in  the  first  15°  as  in  the  reaming  75°. 


a)  Separation  distance  b)  Mach  number 

Figure  56.  Effect  of  diffraction  angle  at  a  point  along  the  bottom  wall  200  mm  downstream  of  the  diffraction 

corner 
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In  cases  D2  and  D4,  the  comer  radius  varied  independently.  The  inerease  in  comer 
radius  from  2.00  mm  to  25.4  mm  had  a  small,  mixed  effect  on  the  separation  distance  and  Maeh 
number  (Fig.  57).  At  a  point  elose  to  the  diffraetion  wall  (x  =25  mm,  y  =  50  mm  relative  to  the 
diffraetion  comer),  the  separation  inereased  with  comer  radius  (Fig.  57a),  and  the  Maeh  number 
decreased  (Fig.  57b).  At  a  point  further  downstream  and  eloser  to  the  bottom  wall  of  the  ehannel 
(x  =  200  mm,  y  =  -50  mm),  separation  deereased,  and  Maeh  number  increased  as  the  comer 
radius  increased.  From  Figure  56,  any  choiee  of  eomer  radius  will  need  to  eonsider  the  loeation 
of  the  refleeting  surfaee  (top  or  bottom  wall)  when  striving  to  optimize  separation  distanee  and 
Maeh  number. 


a)  Separation  distanee  b)  Mach  number 

Figure  57.  Effect  of  corner  radius  at  two  location  (200,  -51)  and  (25,  51) 
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Reflection  angle  and  obstacle  location  results,  cases  R1-R8 

The  results  of  the  cases  D1  through  D4  promote  the  use  of  geometry  from  case  D2  to 
investigate  the  placement  and  angle  of  a  reflecting  surface.  Comer  radius  had  no  significant 
effect  on  the  diffraction,  and  the  90°  diffraction  angle  minimizes  the  transition  length  and  fuel 
requirements.  In  cases  R1-R8,  the  diffraction  angle  was  90°  and  the  comer  radius  was  2  mm 
with  one  exception,  Case  R5.  Case  R5  revisited  the  25.4  mm  comer  radius  to  determine  any 
unforeseen  effect  the  radius  might  have  on  the  combination  of  diffraction  comer  and  reflecting 
surface.  Improvement  is  signified  by  a  reduction  in  separation  and  increase  in  Mach  number  in 
the  R  and  M  series  test  cases.  A  secondary  detonation  will  have  zero  separation  and  a  Mach 
number  in  excess  of  the  CJ  Mach  number  (5.3). 
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Case  R1 


In  case  Rl,  a  45°  reflecting  surface  began  on  the  straight  wall  162.2  mm  downstream  (xq 
=  3.19  hinitiai)  of  the  plane  of  the  diffraction  comer  (Fig.  58a).  Detonations  began  to  diffract  at 
the  comer  and  completely  decoupled  by  the  time  the  shock  reached  the  reflecting  surface  (Fig. 
58b)  as  evidenced  by  the  non-zero  separation  distance.  In  three  of  the  ten  duplicate  mns, 
reinitiation  occurred  at  the  intersection  of  the  ramp  and  the  bottom  wall.  In  every  case  where 
reinitiation  occurred,  the  secondary  diffraction  at  the  end  of  the  ramp  caused  decoupling  of  the 
secondary  detonation.  In  the  cases  where  reinitiation  did  not  occur,  ignition  of  the  reactants 
occurred  and  a  secondary  flame  followed  the  reflected  shock  along  the  surface. 

In  Fig.  58b,  the  low  probability  of  reinitiation  and  quick  decoupling  of  the  resulting 
secondary  detonations  resulted  in  a  small  decrease  in  the  fitted  separation  distance  near  the 
reflecting  surface  compared  to  case  D2.  Likewise,  the  Mach  number  (Fig.  58c)  is  greater  near 
the  reflecting  surface  than  in  case  D2.  The  geometry  in  case  Rl  had  a  30%  chance  to  produce  a 
secondary  detonation  where  all  of  the  D-series  cases  had  a  0%  chance. 
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Figure  58.  Case  R1  (P  =  45°,  xO  =  162  mm,  yO  =  -50.8  mm)  average  of  10  runs 
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Case  R2 


In  case  R2,  a  45°  reflecting  surface  began  on  the  straight  wall  84.7  mm  (1.67  hinitiai) 
downstream  of  the  diffraction  comer  (Fig.  59a).  In  this  arrangement,  the  detonations  did  not 
fully  decouple  before  reaching  the  reflecting  surface  as  evidences  by  the  separation  distance  at 
the  beginning  of  the  reflecting  surface.  Reinitiation  occurred  in  all  of  the  ten  duplicate  mns 
(100%  probability)  and  the  resulting  region  of  small  separation  centered  at  x  =  150  mm,  y  =  25 
mm  is  apparent  in  Figure  59b.  None  of  the  secondary  detonation  waves  survived  the  second 
diffraction  and  the  fit  predicts  decoupling  within  50  mm  of  the  end  of  the  reflecting  surface  (x  = 
135  mm).  The  Mach  number  exceeded  5.0  in  the  region  of  secondary  detonation  (Fig.  59c).  The 
peak,  fitted  Mach  number  was  lower  than  the  CJ  Mach  number  observed  in  Case  D1  (5.3) 
because  of  mn-to-mn  variation  of  the  size  of  the  secondary  detonation.  Case  R2  had  the  highest 
chance  of  reinitiation  with  100%  of  the  duplicate  trials  resulting  in  a  secondary  detonation  and 
the  most  preferable  Mach  number  and  separation  distance  among  the  R-series  cases. 
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Figure  59.  Case  R2  (P  =  45°,  Xq  =  84.7  mm,  yo  =  -50.8  mm)  average  of  10  runs 
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Case  R3 


In  case  R3,  the  placement  of  the  reflecting  surface  was  approximately  the  same  as  in  case 
R2,  but  the  reflecting  angle  was  30°  rather  than  45°  (Fig.  60a).  The  reduced  angle  reduced  the 
compression  of  the  reflected  shock,  and  reduced  the  size  of  the  region  of  zero-separation 
associated  with  reinitiated  detonation  (Fig.  60a).  The  probability  of  reinitiation  was  lower  than 
in  Case  R2  because  seven  of  the  ten  duplicate  runs  successfully  reinitiated  (70%  probability). 
The  decreased  probability  of  reinitiation  increased  the  mean  separation  in  Figure  60b.  As  in 
cases  R1  and  R2,  all  of  the  secondary  detonations  decoupled  fully  after  the  second  diffraction 
comer.  Full  separation  occurs  within  25  mm  of  the  second  diffraction  comer.  Mach  number 
(Fig.  60c)  better  indicates  the  region  of  detonation  above  the  reflecting  surface  though  the  peak 
Mach  number  was  further  reduced  by  the  lower  probability  of  reinitiation  compared  to  case  R2. 
Case  R3  had  the  second  highest  chance  of  reinitiation  among  the  R-series  test  cases.  Only  R2 
had  higher. 
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Figure  60.  Case  R3  (P  =  30°,  Xq  =  80.1  mm,  yo  =  -50.8  mm)  average  of  10  runs. 
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Case  R4 


Case  R4  continued  the  trend  of  decreasing  reflection  angle  (45°  -  30°  -  15°)  (Fig.  61a). 
Due  the  length  of  the  obstacle  needed  for  the  surface  height  to  equal  that  of  cases  R2  and  R3,  the 
reflection  surface  began  directly  opposite  the  diffraction  comer  so  that  it  would  fit  in  the  test 
section.  Unlike  the  higher  reflection  angles,  the  probability  of  reinitiation  was  very  low  with  one 
instance  of  reinitiation  in  six  duplicate  trials  (17  %  probability).  The  fitted  separation  between 
the  shock  and  flame  continued  to  grow  over  time  despite  the  reflection,  and  there  was  no  region 
of  low  separation  (Fig.  61b)  and  high  Mach  number  (Fig.  61c)  typical  of  reinitiation.  The  lack  of 
reinitiation  when  xq  =  0  mm  suggests  that  the  ramp  angle  was  too  shallow  for  reinitiation  at  any 
position  downstream  of  the  diffraction  comer.  With  a  17%  chance  of  reinitiation,  case  R4 
ranked  sixth  of  the  R-series,  better  than  only  R8  and  R7. 
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Figure  61.  Case  R4  (P  =  15°,  Xo  =  0  mm,  yo  =  -50.8  mm)  average  of  6  runs 
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Case  R5 


Case  R5  revisited  the  larger  diffraetion  eomer  radius  (r  =  25.4  mm)  from  ease  D3  (Fig. 
62a).  No  change  relative  to  case  R3  in  the  probability  of  reinitiation,  the  fitted  separation 
distance  (Fig.  62b),  or  the  fitted  Mach  number  (Fig.  62c)  was  expected,  and  the  results  confirm 
that  expectation.  Comparing  case  R3  to  case  R5,  there  was  the  same  small  increase  in  Mach 
number  and  small  reduction  in  separation  distance  near  the  diffraction  corner  as  reported  in  Case 
D3  relative  to  case  D2.  Reinitiation  occurred  in  five  of  the  six  duplicate  trials  (83%  probability), 
which  is  the  same  as  case  three  within  the  margin  of  error  (±17%).  There  appears  to  be  no 
significant  advantage  to  increased  diffraction  comer  radius.  Case  R5  tied  case  R3  for  second 
highest  chance  of  reinitiation  among  the  R-series  cases  and  demonstrates  no  advantage  or 
disadvantage  to  a  larger  diffraction  comer  radius  within  the  range  of  0.2  to  25  mm. 
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Figure  62.  Case  R5  (r  =  25.4  mm,  p  =  30°,  Xo  =  80.1  mm,  yo  =  -50.8  mm)  average  of  5  runs 
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Case  R6 


Beginning  with  case  R6,  the  reflecting  surface  began  one  initial  channel  height  (50.8 
mm)  above  the  diffraction  comer  on  the  same  side  of  the  channel  (Fig.  63a).  In  this 
configuration,  the  Mach  number  before  reflection  was  lower  than  in  cases  R1-R5.  Case  R6  was 
the  first  case  to  exhibit  trapping  of  a  secondary  detonation  between  the  vertical  wall  and  the 
reflecting  surface.  The  beginning  of  the  reflecting  surface  was  the  point  where  diffraction  wall 
and  the  obstacle  met  (xq  =  0).  In  all  of  the  eight  duplicate  mns,  the  decoupled  shock  encountered 
the  obstacle  near  the  end  of  the  reflecting  surface  (xq  =  50  mm),  and  the  separation  distance  was 
too  small  to  allow  a  secondary  detonation  to  propagate  downstream  (Fig.  63b).  In  two  of  the 
rans  reinitiation  occurred  (25%  probability),  but  the  secondary  detonation  waves  were  trapped 
between  the  vertical  wall,  the  reflecting  surface,  and  the  decoupled  combustion  front.  Neither 
the  separation  distance  in  Figure  64b  nor  the  Mach  number  in  Figure  64c  indicated  any 
improvement  (reduced  separation  distance)  over  case  D2  due  to  the  presence  of  the  reflecting 
surface.  Case  R6  demonstrated  a  low  chance  of  reinitiation  (25%)  ranking  fifth  out  of  the  R- 
series  cases,  but  the  geometry  was  completely  ineffective  because  the  secondary  detonations 
formed  in  a  region  of  reactants  bounded  by  the  reflecting  surface  and  the  decoupled  combustion 
front.  The  infrequent,  secondary  detonations  quenched  once  the  reactants  in  that  region  were 
consumed. 
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Figure  63.  Case  R6  (P  =  135°,  Xo  =  0  mm,  yo  =  50.8  mm)  average  of  8  runs. 
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Case  R7 


In  case  R7,  the  reflecting  surface  was  moved  to  xq  =  43  mm  (0.85  hinitiai)  on  the  upper 
wall  of  the  channel  (Fig.  64a).  At  this  location,  the  first  contact  between  the  decoupled  shock 
and  the  reflecting  surface  occurred  near  the  midpoint  of  the  reflecting  surface.  The  incident 
Mach  number  at  the  reflecting  surface  was  the  lowest  of  any  case  yet  and  the  resulting  shock 
reflection  was  insufficient  for  reinitiation.  Reinitiation  occurred  in  none  of  the  eight  duplicate 
runs  (0%  probability).  The  separation  distance  (Fig  64b)  and  Mach  number  (Fig.  64c)  exhibit  no 
significant  difference  from  case  R6  or  from  case  D2.  As  in  case  R6,  a  pocket  of  reactants 
became  isolated  between  the  top  wall  of  the  channel,  the  reflecting  surface  and  the  combustion 
front,  but  unlike  case  R6  there  was  no  secondary  detonation  to  trap  in  case  R7.  Case  R7 
performed  the  worst  of  the  R-series  cases  with  no  chance  of  reinitiation  due  to  low  shock  Mach 
number  at  the  reflecting  surface. 
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Figure  64.  Case  R7  (P  =  135°,  Xo  =  43  mm,  yo  =  50.8  mm)  average  of  8  runs. 
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Case  R8 


In  the  final  reflection  case,  the  reflecting  surface  was  located  at  xq  =  169  mm  (3.33  hinitiai) 
on  the  top  wall  of  the  channel  (Fig.  65a).  The  main  difference  between  this  and  case  R7  was  the 
reflection  of  the  decoupled  shock  wave  from  the  upper  wall  of  the  channel  before  encountering 
the  reflecting  surface.  The  combination  of  two  reflections,  one  from  the  top  of  the  channel  and 
one  from  the  reflecting  surface,  was  sufficient  for  reinitiation  to  occur  in  four  of  the  nine 
duplicate  runs  (44%  probability).  Along  the  reflecting  surface,  the  separation  distance  decreased 
(Fig.  65b),  and  the  Mach  number  increased  (Fig.  65c)  compared  to  case  D2.  In  every  case,  the 
secondary  detonation  decoupled  after  the  second  diffraction  comer.  How  far  the  secondary 
detonation  traveled  before  fully  decoupling  is  unknown  because  the  secondary  detonations  were 
still  partially  coupled  when  they  reached  the  end  of  the  image  frame.  This  configuration  had  the 
best  probability  of  reinitiation  among  the  geometries  with  a  positive  vertical  offset  (yo  =  hinitiai), 
but  there  was  no  configuration  tested  wherein  reinitiation  occurred  100%  of  the  time.  Case  R8 
had  the  second  lowest  change  of  reinitiation  of  the  R-series  cases  and  performed  better  than  case 
R7  only  because  of  the  additional  shock  reflection  from  the  top  of  the  channel.  Based  on  cases 
R6-R8,  vertical  offset  should  be  minimized  in  the  design  of  a  detonation  diffuser. 
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Figure  65.  Case  R8  (P  =  135°,  Xo  =  169  mm,  yo  =  50.8  mm)  average  of  9  runs. 
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The  cases  in  the  R-series  explored  a  parameter  space  consisting  of  reflecting  angle  and 
location  to  determine  where  reinitiation  was  likely  to  occur.  Figure  66  shows  how  the 
probability  of  reinitiation  depends  on  the  reflecting  surface  angle  and  location.  Reinitiation  was 
certain  in  case  R2,  and  reinitiation  occurred  intermittently  in  cases  Rl,  R3-R5,  and  R8.  Figure 
66a  shows  that  the  probability  of  reinitiation  increases  with  reflection  angle  (P).  As  P  increases, 
the  component  of  the  shock  velocity  normal  to  the  surface  increases  resulting  in  higher 
compression  and  higher  probability  of  reinitiation.  The  probability  of  reinitiation  as  a  function 
of  xo  in  Fig.  66b  depended  greatly  on  the  vertical  offset,  yo.  As  the  vertical  offset  increased  from 
-hinitiai  to  hinitiai,  the  minimum  xq  for  reinitiation  increased  quickly,  beginning  at  0.0  and 
increasing  to  1.67  hinitiai  at  a  vertical  offset  of  yo  =  hinitiai-  The  minimum  limit  of  xq  for 
reinitiation  occurred  because  secondary  detonations  were  trapped  by  the  decoupled  combustion 
fronts.  As  xq  increased,  larger  separation  distances  allowed  secondary  detonations  to  escape  the 
combustion  fronts,  but  a  single  shock  reflection  was  insufficient  to  reinitiate  detonation.  Two 
reflections  of  the  decoupled  shock  (once  from  the  top  wall  of  the  channel  and  once  from  the 
reflecting  surface)  were  sufficient  to  increase  the  probability  of  reinitiation  to  17%  at  xq  = 
3.33  hinitiai-  Once  xq  was  sufficiently  large  (xq  >  1-67  hinitiai),  reinitiation  was  no  longer  certain 
(100%  probability)  on  the  bottom  of  the  channel  (yo  =  -hinitiai)  because  the  decoupled  shock 
decayed  too  much  before  encountering  the  reflecting  surface- 
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Figure  66.  Trends  in  the  probability  of  reinitiation 


98 


Iterative  multiple  obstacle  case  results,  cases  Ml-Mll 

The  final  set  of  test  cases  expands  upon  the  quantitative  data  from  the  reflection  cases. 
The  number  of  different  test  cases  in  the  multiple  obstacle  set  made  quantitative  analysis 
impractical,  and  it  was  unnecessary  for  any  case  without  successful  transmission  of  a  secondary 
detonation.  Of  the  twelve  cases,  none  resulted  in  detonation  in  the  supercritical  channel,  but 
each  contributed  to  understanding  of  the  secondary  detonation  and  the  second  diffraction  comer. 
The  second  diffraction  comer  remains  the  obstacle  to  a  successful  detonation  diffuser. 

Case  Ml 

In  case  Ml,  the  second  diffraction  comer  had  a  radius  of  13  mm  in  an  attempt  to  prevent 
decoupling  (Fid.  67a).  Since  the  data  analyses  from  the  two  previous  cases  with  large  diffraction 
comer  radii  were  incomplete  at  the  time,  this  was  considered  a  useful  test  case.  In  each  of  four 
trials,  diffraction  occurred  after  the  obstacle  and  the  secondary  detonation  failed  (Fig.  67b). 
Figure  66b  is  a  composite  of  frames  from  a  single  video  highlighting  the  reinitiation  and 
subsequent  diffraction.  The  shock  and  flame  in  each  frame  was  shaded  in  a  different  color  to 
distinguish  that  frame  from  its  neighbors.  The  chronological  order  of  shading  colors  is  blue, 
green,  red,  and  black.  After  four  frames,  the  cycle  repeats. 
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a)  Schematic 
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Figure  67.  Case  Ml,  13  mm  high  obstacle  with  rounded  diffraction  corner 
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=  241  mm 


Case  M2 


In  case  two,  reinitiation  began  in  two  places  (Fig.  68b).  The  first  was  either  the  third  or 
fourth  obstacle  and  the  second  was  one  of  the  last  three  obstacles.  Neither  secondary  detonation 
traveled  across  the  shock  front  without  decoupling.  After  this  case,  the  sharp  diffraction  angle  at 
the  top  of  each  ramp  was  rounded  to  reduce  the  severity  of  the  diffraction.  The  reduced  height  of 
the  obstacles  continued  to  have  no  effect  on  reinitiation. 


a)  Schematic 


b)  Composite  Frame 
At  =  20.4  ps 


Figure  68.  Case  M2:  multiple  6.4  mm  high  ohstacles 
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Case  MS 


The  rounded  diffraetion  comers  increased  the  number  of  reinitiation  events  in  case  three 
from  two  to  three  (Fig.  69b).  Again,  the  secondary  detonations  decoupled  before  reaching  the 
upper  wall  of  the  channel.  After  case  M3,  the  channel  height  was  reduced  to  102  mm  to  decrease 
the  distance  the  secondary  detonation  waves  and  their  associated  decoupled  shocks  had  to  travel 
before  encountering  a  solid  surface. 


a)  Schematic 


b)  Composite  Frame: 
At  =  20.4  ps 


Figure  69.  Case  M3:  multiple  6.4  mm  high  obstacles  with  rounded  corners 
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Case  M4 


The  smooth  upper  wall  of  the  shortened  channel  in  Case  M4  (Fig.  70a)  reflected  the 
decoupled  shocks,  but  the  reflections  were  too  weak  for  reinitiation  (Fig.  70b).  The  transverse 
shocks  from  the  second  and  third  secondary  detonations  did  not  reach  the  upper  wall  before  the 
leading  shock  reached  the  end  of  the  test  section.  In  the  next  iteration,  a  second  set  of  obstacles 
was  added  to  the  top  of  the  channel  to  form  more  transverse  waves  in  order  to  foster  more 
secondary  detonations. 


Figure  70.  Case  M4:  multiple  6.4  mm  high  obstacles  in  102  mm  tall  channel 
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Case  M5 


With  two  sets  of  obstacles  (Fig.  71a),  the  picture  became  muddled  with  several  sets  of 
transverse  waves  overlapping  the  lead  shock  in  previous  frames  (Fig.  71b).  Two  or  three 
secondary  detonations  formed  on  the  bottom  of  the  channel,  but  none  formed  on  the  top.  This 
was  due  to  the  low  Mach  number  of  the  lead  shock  at  the  top  of  the  channel.  The  separation 
distance  in  the  final  frame  was  the  smallest  of  the  cases  so  far.  In  the  run  shown  in  figure  71b 
part  of  the  wave  front  remained  coupled,  but  that  was  not  consistent  across  multiple  duplicate 
trials. 
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Figure  71.  Case  M5:  102  mm  high  channel  with  obstacles  on  both  walls 


104 


Case  M6 


To  better  foster  secondary  detonation  at  the  top  of  the  channel,  the  obstacles  at  the  top 
had  their  diffraction  radius  increased  in  Case  M6  (Fig.  72a).  Partially  coupled  wave  fronts  at  the 
end  of  the  test  section  were  rare  and  no  more  likely  than  in  Case  M4.  Multiple  transverse  shocks 
continued  to  muddle  the  composite  though  there  were  no  more  secondary  detonations.  There 
was  no  significant  difference  in  the  behavior  in  the  two  cases.  At  this  point,  the  test  cases 
changes  to  a  different  avenue,  modifying  the  diffraction  comer  in  conjunction  with  a  series  of 
reflection  obstacles. 


a)  Schematic 
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Figure  72.  Case  M6:  rounded  obstacles  on  top  and  bottom  walls 
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Case  M7 


In  case  M7  (Fig.  73a),  the  stepped  diffraction  comer  caused  one  secondary  detonation, 
and  the  obstacles  on  the  top  of  the  channel  caused  none  in  any  of  the  repeat  trials.  This 
configuration  performed  poorly  compared  to  the  case  M6,  and  the  separation  distance  was 
clearly  larger  in  the  final  frame  (Fig.  73b).  This  case  was  the  first  to  exclude  obstacles  on  the 
bottom  of  the  channel.  In  the  next  case,  the  obstacle  set  was  placed  on  the  bottom  of  the  channel. 


a)  Schematic 
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b)  Composite  Frame: 
At  =  20.4  ps 


Figure  73.  Case  M7:  Stepped  diffraction  with  obstacles  on  top  wall 
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Case  M8 


Moving  the  set  of  reflecting  obstacles  to  the  bottom  of  the  channel  was  an  improvement 
(Fig.  74a)  because  reinitiation  occurred  at  least  once.  There  was  at  least  one  reinitiation  on  the 
stepped  diffraction  (Fig.  74b),  and  none  among  the  obstacle  set.  At  this  point,  it  was  clear  that 
the  stepped  diffraction  comer  performed  worse  than  previous  cases  and  the  next  case  increased 
the  size  of  the  steps  and  put  some  space  between  the  obstacles  in  order  to  create  two  reflections 
per  diffraction. 


a)  Schematic 
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Figure  74.  Case  M8:  Stepped  diffraction  with  obstacles  on  bottom  wall 
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Case  M9 


Because  decoupling  occurred  despite  the  stepped  diffraction  comer,  it  seemed  possible  to 
reinitiate  detonation  using  the  decoupled  shock.  In  case  M9,  the  steps  were  redesigned  to 
promote  reinitiation  by  splitting  the  single  90°  turn  in  the  wall  into  two  45°  reflections  (Fig.  75). 
The  size  of  the  steps  was  increased  to  reduce  the  number  of  diffractions  that  a  secondary 
detonation  would  encounter,  and  a  similar  set  of  obstacles  was  placed  on  the  bottom  of  the 
channel  (Fig.  76a).  In  this  configuration,  reinitiation  always  occurred  at  the  first  and  second 
obstacles  on  the  bottom  of  the  channel  but  never  at  the  third  obstacle  (Fig.  76b).  The  first 
obstacle  along  the  diffracting  wall  caused  reinitiation  in  every  case,  but  the  second  never  did. 
The  first  bottom  obstacle  constricted  the  channel  height  in  such  a  way  that  the  combined 
diffraction  of  the  initial  diffraction  comer  and  the  diffraction  comer  at  the  top  of  the  obstacle 
reduced  the  shock  strength  too  much  for  any  further  secondary  detonations.  This  was 
unexpected  at  the  second  obstacle  on  the  bottom  of  the  channel  should  also  cause  a  reinitiation 
based  on  reflection  case  two.  The  obstacle  also  restricted  flow  during  filling  of  the  test  section 
forcing  increased  fill  pressure  an  undesired  side  effect. 
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Figure  75.  Two-reflection  geometry  for  diffraction  step 
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a)  Schematic 


b)  Composite  Frame: 
At  =  48.8  ps 


Figure76.  Case  M9:  Double  reflection  with  restriction 
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Case  MIO 


To  relieve  the  flow  restriction  caused  by  the  first  obstacle,  and  stagger  the  diffraction  of 
the  initial  comer  and  the  obstacles,  the  first  obstacle  was  removed  for  case  MIO  (Fig.  77a).  The 
first,  bottom  obstacle  caused  secondary  detonation  in  all  five  duplicate  trials  (Fig.  77b).  The 
secondary  wave  decoupled  before  traversing  the  channel  in  each  case,  and  none  of  the  other 
obstacles  cause  reinitiation.  Unlike  case  M9,  some  of  the  trials  saw  shock  ignition  on  the  top 
two  obstacles  and  on  the  second  bottom  obstacle.  Elevated  fill  pressures  were  unnecessary  with 
the  first  obstacle  removed. 


no 


a)  Schematic 


b)  Composite  Frame: 
At  =  48.8  ps 


Figure  77.  Case  MIO:  Double  reflection  shape  without  restriction 
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Case  Mil 


The  final  case  arranged  four  obstacles  so  that  the  lead  shock  encountered  all  of  them 
within  20|as  (Fig78).  The  obstacles  were  placed  with  the  leading  edges  on  a  circle  of  radius  50.8 
mm  centered  opposite  the  initial  diffraction  corner  (Fig  79a).  The  reflection  angle  of  each 
obstacle  exceeded  45°  relative  to  the  incident,  decoupled  shock  wave,  and  the  incident  shock 
Mach  number  exceeded  necessary  for  reinitiation  obtained  from  the  R-series  test  cases.  Case 
Mil  was  the  first  configuration  to  utilize  the  shock  Mach  number  data  from  case  D2  and  the 
reinitiation  requirements  obtained  in  the  R-series  test  cases.  This  critical  information  was  not 
available  earlier  because  of  the  long  data  processing  times  required  obtain  the  shock  Mach 
number  in  case  D2. 

The  result  of  the  informed  design  of  case  Mil  was  four  secondary  detonation  waves,  and 
no  decoupling  of  the  initial  detonation  (Fig.  79b).  The  five  detonation  waves  traveled  down  the 
expanding  channels  (0i  =  10°)  between  the  obstacles.  Often,  partial  decoupling  occurred,  but  all 
of  the  waves  reached  the  ends  of  the  obstacles  at  about  the  same  time.  The  secondary  diffraction 
comers  completely  decoupled  the  secondary  detonations,  but  the  collision  of  decoupled  shocks 
leaving  two  adjacent  channels  sometimes  caused  another  secondary  detonation.  This  was  the 
best  result  of  the  multiple  obstacle  cases  because  it  exhibited  the  spontaneous  reinitiation  that 
defines  detonation  diffraction  from  critical  channels. 
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Figure  78  Image  sequence  of  case  Mil:  Run  5  At  =  20  jis. 
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a)  Schematic 


b)  Composite  Frame: 
At  =  48.8  ps 


Figure  79.  Case  Mil:  Split  channel  geometry 
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Overall  trends  in  reinitiation  for  the  M  series  cases 

The  multi-obstacle  cases  utilized  different  numbers  of  obstacles,  and  in  general,  the  cases 
with  more  obstacles  had  more  reinitiations.  The  added  obstacles  come  at  the  cost  of  higher  flow 
loss.  Figure  78  shows  that  each  additional  obstacle  has  a  reduced  maximum  chance  of  causing  a 
reinitiation.  For  14  or  more  obstacles,  the  chance  of  reinitiation  is  at  most  86%.  Based  on 
Figure  79,  no  more  than  12  obstacles  should  be  used  for  the  current  channel  height  and  fuel 
because  there  is  no  additional  benefit  from  the  added  obstacles. 


Figure  80.  Diminishing  return  of  additional  obstacles 
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For  a  simple  quantitative  comparison  of  the  different  configurations,  the  separation 
distance  was  measured  at  two  points  in  each  case.  One  point  was  at  the  bottom  of  the  channel  in 
the  last  frame  where  the  leading  shock  was  visible,  and  the  other  was  along  the  diffraction  wall 
in  the  frame  before  the  lead  shock  encountered  the  top  of  the  channel.  These  measurements 
provided  a  basic  indicator  of  the  performance  in  each  case.  Smaller  separations  were  preferred 
since  they  indicated  the  least  time  after  decoupling.  Figure  80  shows  how  the  separation 
distances  evolved  through  the  series.  In  the  early  cases  (Ml -M3),  low  separation  at  the  bottom 
of  the  channel  came  with  high  separation  at  the.  Starting  with  case  M4,  separation  at  the  top 
began  to  match  the  bottom,  and  at  the  end  of  the  series,  both  are  at  their  lowest. 

50 


Case  Number 


■  Bottom  Wall  ■  Diffraction  Wall 


Figure  81.  Comparison  of  final  separation  distance  for  M-series  cases 
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VI.  Discussion 


Crossover  tube  studies  identify  diffraction  and  reflection  parameters. 

In  a  crossover  tube,  the  incoming  detonation  first  diffracts  when  it  encounters  the 
entrance  of  the  crossover  tube  (Fig.  15).  Then  the  decoupled  shock  reflects  from  the  reflecting 
surface  reinitiating  detonation.  Shock  reflection  and  subsequent  reinitiation  of  detonation 
inspired  the  concept  of  a  detonation  diffuser  that  used  a  combination  of  diffraction  and 
reinitiation  to  transmit  a  detonation  from  a  subcritical  channel  to  a  supercritical  channel  without 
going  through  DDT. 

Crossover  tube  experiments  were  useful  for  identifying  some  of  the  important  parameters 
in  a  detonation  diffuser  (Figs.  16  and  17).  The  diffraction  angle,  diffraction  comer  radius,  and 
the  reflecting  surface  angle  were  the  three  parameters  that  influenced  the  diffraction  and 
reinitiation  in  crossover  tubes.  Each  parameter  was  studied  in  depth  to  determine  the  effect  on 
reinitiation. 

Diverging/Converging  experiment  establishes  feasibility  and  benefit. 

To  determine  whether  a  detonation  diffuser  utilizing  decoupling  and  reinitiation  was 
feasible  and  beneficial,  an  experiment  was  carried  out  to  compare  a  diverging  geometry  to  a  step 
expansion  followed  by  a  converging  wall  (Figs.  50  and  51).  The  divergence  angle  in  the 
diverging  case  was  less  than  the  theoretical  maximum  for  a  detonation  to  remain  coupled,  but  in 
disagreement  with  theory  decoupling  was  observed  in  schlieren  video.  In  the  converging  case, 
the  primary  detonation  decoupled  before  it  encountered  the  converging  wall.  Then  a  large 
increase  in  chemiluminescence  near  the  midpoint  of  the  converging  section  signaled  reinitiation. 
Analysis  of  the  two  cases  indicated  that  the  converging  configuration  would  result  in  a  2.3% 
shorter  transition  to  a  supercritical  channel  height  than  a  diverging  channel. 
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Diffraction  cases  indicate  small  diffraction  angle  preferred  and  corner  radius  trends  mixed. 

Following  the  Diverging/Converging  experiment,  the  results  aim  to  satisfy  the  research 
objective  defined  in  Chapter  1.  The  diffraction  test  cases  (D-series)  provide  the  required  maps  of 
shock  Mach  number  and  separation  distance  (Figs.  53-55).  In  the  D-series  of  test  cases, 
measurements  of  the  separation  distance  and  shock  Mach  number  indicated  that  detonations 
partly  decouple  at  a  diffraction  angle  of  15°  and  fully  decouple  at  larger  angles  (Figs.  53  and  54). 
Partial,  temporary  decoupling  explains  the  decoupling  that  was  observed  at  diffraction  angles 
less  than  the  theoretical  limit  in  the  previously  discussed  diverging/converging  experiment.  The 
separation  distance  increases  linearly  with  diffraction  angle  at  a  point  200  mm  (3.94  hinitiai) 
downstream  of  the  diffraction  comer  (Fig.  56a).  The  Mach  number  also  decreased  with 
increasing  diffraction  angle,  but  the  trend  was  nonlinear  with  larger  rate  of  decline  between  0° 
and  15°  than  between  15°  and  90°  (Fig.  56b).  Increasing  the  comer  radius  had  mixed  results 
depending  on  location.  Near  the  bottom  of  the  channel,  the  separation  distance  decreased  and 
the  Mach  number  increased  (both  of  which  are  preferred)  (Fig.  57a).  At  a  point  along  the 
vertical  wall  (yo  =  hinitiai),  the  separation  increased  and  the  Mach  number  decreased  at  rates 
similar  to  the  improvements  along  the  bottom  of  the  channel  (Fig.  57).  Later,  in  case  R5  the 
comer  radius  was  shown  to  have  no  statistically  significant  effect  on  reinitiation  (Fig.  62).  The 
combination  of  mixed  and  insignificant  effects  suggests  that  an  arbitrary  choice  of  comer  radius 
is  acceptable  in  the  range  of  2-25  mm. 

Chance  of  reinitiation  depends  strongly  on  reflection  angle  and  position. 

The  R-series  test  cases  found  the  range  in  each  of  three  parameters  (P,  xq,  and  yo)  were 
reinitiation  occurred  and  the  probability  of  reinitiation  within  those  ranges.  Reinitiation  of 
detonation  via  oblique  shock  reflection  had  two  operational  regions.  In  the  first  region. 
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reinitiation  was  certain  due  to  ample  compression  of  reactants  and  sufficient  separation  between 
the  shock  and  combustion  fronts  (Fig.  66a).  Reinitiation  was  certain  for  the  operating  space 
where  |3  =  45°,  yo  =  -hinitiai,  and  xO  <  1.67  hinitiai  (Fig.  66a).  In  the  second  region,  reinitiation 
occurred  in  some  fraction  of  the  repeated  trials  (Fig.  66a).  The  probability  of  reinitiation 
increased  with  increasing  |3,  decreased  with  increasing  xq  and  decreased  with  increasing  yo  (Fig. 
66).  When  yo  =  hinitiai,  reinitiation  only  occurred  for  xo  =  3.33  hinitiai  and  only  44%  of  the  time 
(Fig.  66b).  There  was  a  tradeoff  between  shock  strength  and  separation  distance  when  yo  =  hintiai- 
At  lower  xo  secondary  detonations  became  trapped  by  the  decoupled  combustion  front,  and  at 
high  Xo  the  lead  shock  decayed  too  much  for  a  secondary  detonation  to  form  (Fig. 66b). 
Multi-obstacle  cases  bridged  the  gap  between  subcritical  and  critical  diffraction  behavior 

The  multi-obstacle  test  cases  satisfied  the  phase  3  research  objectives  by  quickly  iterating 
on  the  size,  number  and  position  of  obstacles  using  a  qualitative  evaluation  of  each  geometry. 
The  M-series  cases  utilized  the  information  gained  from  the  D  and  R  series  to  bridge  the  gap 
between  subcritical  and  critical  diffractions.  Two  cases  (M6  and  Mil)  maintained  partial 
coupling  when  the  wave  reached  the  end  of  the  test  section  (Figs.  72  and  78).  One  of  the  two 
(Mil)  also  exhibited  spontaneous  reinitiation  after  the  obstacles  (Fig.  78b).  In  case  Mil,  the 
spontaneous  reinitiation  of  detonation  indicated  critical  diffraction  behavior.  The  results  indicate 
that  a  fully  successful  detonation  diffuser  should  initiate  multiple  secondary  detonations  that 
combine  to  form  a  fiilly  coupled  detonation  front.  The  first  two  steps  for  improving  the 
geometry  in  following  design  iterations  are  to  eliminate  the  expansion  in  the  channels  between 
obstacles  and  to  reduce  the  diffraction  angle  at  the  end  of  the  obstacles. 
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VII.  Conclusion  and  Future  Work 


A  detonation  wave  in  a  subcritical  channel  decouples  when  diffracted.  Decoupling  is 
undesired  in  the  transition  between  a  predetonator  and  the  thrust  tube  in  PDEs  where  the  smallest 
possible  predetonator  minimizes  weight  and  fuel  requirements.  A  detonation  diffuser  of  the  type 
studied  in  this  research  reinitiates  detonation  after  the  decoupling  by  reflecting  the  decoupled 
shock  wave  back  into  itself.  The  initial  diffraction  angle,  0,  should  be  as  small  as  space 
considerations  allow,  and  the  comer  radius,  r,  can  be  anywhere  between  2  and  25  mm  with  no 
adverse  effect.  An  oblique  shock  reflection  is  sufficient  to  reinitiate  detonation  with  sufficient 
shock  strength  and  mixture  sensitivity.  The  reflection  angle  must  be  45°  or  greater  for  certain 
reinitiation.  The  lead  obstacle  if  placed  on  the  wall  of  the  channel  opposite  the  diffraction  comer 
(bottom  wall)  should  be  no  more  than  1.67  hinitiai  downstream.  Moving  the  lead  obstacle  away 
from  the  bottom  wall  reduced  the  minimum  downstream  distance  and  the  minimum  distance  was 
0  when  the  lead  obstacle  was  hinitiai  above  the  diffraction  comer.  Multiple  obstacles  in  series 
caused  multiple  reinitiations,  but  there  was  a  diminishing  return  associated  with  each  additional 
obstacle.  No  more  than  12  obstacles  should  be  used  in  series.  A  better  design  reinitiated 
detonation  in  four  separate  sub  channels,  and  exhibited  the  spontaneous  reinitiation  that  defines  a 
critical  diffraction. 

Future  research  in  detonation  diffraction  and  reinitiation  should  concentrate  on  separate 
reinitiation  of  several  secondary  detonations  from  the  initially  decoupled  shock.  The  first  two 
steps  of  further  investigation  are  to  reduce  the  diffraction  angle  at  the  ends  of  the  obstacles  in 
case  Mil  and  to  eliminate  the  expansion  in  the  channels  between  obstacles.  The  ultimate  goal  is 
to  achieve  fully  coupled  detonation  at  a  super-critical  channel  height  at  the  exit  of  the  diffuser. 
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Once  the  channel  height  exceeds  the  critical  height,  no  special  geometry  is  necessary  to  expand 
into  arbitrarily  large  thrust  tubes. 
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IX.  Appendix  A  -  Schlieren  technique,  equipment,  and  uncertainty 
Technique 

Schlieren  visualization  takes  advantage  of  the  coupling  between  density  and  refractive 
index  to  make  structures  visible  in  transparent  media  (Settles,  2001). 

condenser 


Figure  A-1.  Z-type  schlieren  arrangement  (Settles,  2001) 

In  Fig.  A-1,  the  primary  mirror  collimates  light  from  a  point  source  (the  slit).  The 
parallel  light  passes  through  the  test  region  where  density  gradients  diffract  some  of  the  rays. 
The  secondary  mirror  focuses  the  remaining  parallel  beam  to  an  image  of  the  source.  The 
diffracted  rays  focus  to  points  away  from  the  image  of  the  source.  The  direction  a  ray  is 
displaced  depends  on  the  sign  of  the  density  gradient,  and  the  distance  depends  on  the  magnitude 
of  the  gradient.  Beyond  the  source  image,  an  image  of  the  test  region  forms.  Without  further 
interference,  structures  in  the  flow  are  visible  as  shadows  with  adjacent  bright  areas  on  either 
side  (Fig.  A-2a).  Adding  a  knife-edge  at  the  source  image  blocks  some  of  the  rays  eliminating 
the  bright  band  on  one  side  of  a  structure  and  revealing  the  direction  of  gradients. 
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Figure  A-2.  Series  of  schlieren  photos  of  a  turbulent  gas  jet  with  increasing  cut-off.  The  cut-off  degree  is  a) 
0%,  b)  20%,  c)  40%,  d)  60%,  e)  80%,  0  90%,  g)  95%,  and  h)  100%.  Photos  by  Rosanna  Quinones  (Settles, 

2001) 


The  selection  of  a  camera  determines  the  spatial  and  temporal  characteristics  of  the 
recorded  images.  The  schlieren  technique  is  analog  and  the  magnification,  resolution,  exposure, 
and  frame  rate  of  the  camera  determine  the  sampling  and  uncertainty  of  the  image. 
Magnification  and  pixel  count  set  the  spatial  resolution.  Frame  rate  and  exposure  set  the 
temporal  resolution. 

Equipment 

The  schlieren  system  used  in  this  research  is  an  adaptation  of  the  z-type  configuration 
(Fig.  A-3).  The  layout  changed  significantly  after  Nielsen’s  crossover  study  (Fig.  A-4)  because 
the  camera  was  too  far  from  the  test  section  to  bring  objects  into  sharp  focus,  and  to  decrease 
distortion  of  the  parallel  beam  due  to  the  large  turning  angles  at  the  fold  mirrors.  The  drawback 
of  the  current  layout  is  that  all  four  mirrors,  the  light  source,  and  the  camera  table  must  be  moved 
to  image  a  different  section  of  the  thrust  tube. 
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Figure  A-3.  Schlieren  arrangement  for  detonation  diffuser  study 


Figure  A-4.  Previous  arrangement 
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The  system  is  eomposed  of  a  eustom  light  souree,  two  318  mm  diameter,  foeusing 
mirrors,  two  318  mm  diameter  flat  mirrors,  a  knife  edge,  and  a  camera. 

The  light  source  consists  of  a  90001m  LED,  two  lenses,  a  pinhole,  and  a  blackout  tube  all 
mounted  to  an  optical  breadboard  (Fig.  A-5).  The  lenses  condense  the  light  from  the  aperture 
diameter  of  the  lamp  to  the  diameter  of  the  adjustable  slit.  The  slit  was  useful  for  alignments 
since  it  produces  a  small  image  at  the  focal  plane,  but  no  aperture  was  used  for  imaging.  The 
pinhole  sharpens  the  focus.  A  blackout  tube  blocks  stray  light  preventing  interference  with  the 
parallel  beam.  The  entire  apparatus  attaches  to  a  breadboard  to  preserve  the  alignment  between 
uses. 


adjustable  150  mm  lens  tube  50  mm 


Figure  A-5.  Light  source  assembly 
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The  focusing  mirrors  are  aluminum  first  surface  mirrors  with  2.54  m  local  length  (Fig.  A- 
6).  The  mirror  substrate  is  Pyrex  ®  ground  to  a  parabolic  curvature  with  focal  length  accuracy  of 
±1.5%.  The  surface  accuracy  is  l/8th  of  the  median  wavelength.  For  the  visible  spectrum,  the 
accuracy  is  70  nm.  The  reflective  coating  is  a  thin  layer  of  pure  aluminum.  A  275  nm  thick  SiO 
layer  prevents  oxidation  and  protects  the  aluminum.  The  adjustable  mounts  hold  the  mirrors  and 
sheet  metal  covers  prevent  damage  between  uses.  The  mounts  bolt  to  heavy  stands  with 
adjustable  legs  for  leveling.  The  flat  mirrors  are  also  high  quality  first  surface  mirrors.  The 
construction  of  the  flat  mirrors  is  identical  to  the  focusing  mirrors.  They  also  have  l/8th 
wavelength  surface  accuracy.  The  mounts  and  stands  used  for  the  flat  mirrors  are  the  same  as  for 
the  focusing  mirrors. 


Figure  A-6.  Focusing  mirror  with  cover 

The  knife-edge  is  an  ordinary  razor  blade  mounted  on  a  translation  stage.  An  optical 
filter  mount  that  bolts  to  the  translation  stage  holds  the  razor  blade.  A  filter  mount  holds  the 
blade  vertically  or  horizontally  to  image  horizontal  or  vertical  gradients  respectively.  The  stage 
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allows  25  mm  of  overall  travel,  and  has  a  fine  threaded  screw  for  fine  placement.  Daily 
inspection  of  the  edge  ensured  that  there  were  no  nicks  or  corrosion  that  could  affect  the  cut-off. 

The  camera  used  for  detonation  diffuser  study  was  a  Phantom  v710.  The  V710  has  a 
1280  by  800  pixel  array  and  full  resolution  frame  rate  of  7530  frame/s.  At  reduced  resolution, 
the  maximum  frame  rate  increases  1.4  million  times  per  second.  The  exposure  is  adjustable  from 
296  ns  to  1  ms.  An  external  trigger  starts  recording  within  100  ms  of  the  first  detonation  in  a 
burst.  A  laptop  equipped  with  control  software  adjusts  the  camera  settings  and  saves  images 
over  a  network  connection. 

Uncertainty 

The  selection  of  a  camera  determines  the  spatial  and  temporal  uncertainty  in  schlieren 
images.  The  schlieren  technique  is  analog,  and  the  camera  governs  magnification,  pixel  count, 
exposure,  and  frame  rate.  Magnification  and  pixel  count  set  the  spatial  resolution  while  frame 
rate  and  exposure  set  the  temporal  resolution. 

As  an  example,  consider  detonation  of  stoichiometric  hydrogen/air  at  standard 
temperature  and  pressure.  Hydrogen  is  widely  used  in  detonation  study  due  to  high  sensitivity 
and  small  cell  size.  High  sensitivity  makes  detonation  easy  to  achieve,  and  small  cell  size 
reduces  the  size  of  experiments.  Stoichiometric  hydrogen/air  has  a  theoretical  CJ  velocity  of 
1971  m/s  (Schultz,  1999),  and  cell  size  of  8.19mm  (Ciccarelli  et  al.,  1994).  In  the  proposed 
optical  test  section,  the  initial  channel  height  is  50.8  mm  or  6.20  X  significantly  less  than  the 
critical  channel  height.  Assuming  a  cell  of  equal  width  and  length,  triple  point  collisions  take 
place  every  4.16  ps.  In  order  to  satisfy  the  Nyquist  Sampling  Theorem  for  cell  size,  the  temporal 
resolution  must  be  shorter  than  2.08  ps,  and  the  spatial  resolution  must  be  smaller  than  4.09  mm. 
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Converting  to  frame  rate,  the  minimum  allowable  frame  rate  is  481000  frame/s.  The  minimum 
number  of  pixels  needed  to  resolve  a  cell  in  the  203  mm  wide  viewing  section  is  50. 

The  frame  rate  and  pixel  count  requirements  vary  depending  on  cell  size  and  CJ  velocity. 
The  minimum  cell  size  for  subcritical  diffraction  in  the  optical  section  is  5.08mm  (H  =  10  X),  and 
the  maximum  is  50.8  mm  (H  =  X).  The  corresponding  spatial  resolutions  are  80  pixels  and  4 
pixels  respectively.  CJ  velocities  range  from  1100  m/s  to  2200  m/s  for  hydrogen/air  depending 
on  equivalence  ratio.  The  corresponding  minimum  frame  rates  are  43000  frame/s  and  874000 
frame/s  respectively.  The  stoichiometric  example  falls  roughly  half  way  between  the  limits.  The 
Phantom  v7.0  camera  used  in  Stevens  et  al.  (2011)  ran  at  265  pixels  by  64-pixel  resolution  and 
75000  frame/s.  As  a  result,  the  temporal  resolution  was  insufficient  to  capture  the  cellular 
structure. 

Four  commercially  available  cameras  meet  the  frame  rate  requirement  for  stoichiometric 
hydrogen/air  detonation.  They  are  the  Phantom  vl2.1  and  Phantom  v710  by  vision  research,  the 
Fastcam  SA-5  by  Photronics,  and  the  HPV-2  by  Shimadzu.  The  HPV-2  has  memory  for  only 
100  frames  making  it  inappropriate  for  observing  the  entire  process  of  diffraction  and  reinitiation 
in  a  detonation  diffuser.  The  remaining  cameras  force  a  trade-off  between  frame  rate  and  image 
resolution.  Table  1  compares  the  resolution  and  frame  rate  settings  of  each  camera.  The 
Phantom  v7.1  is  included  for  reference.  The  Phantom  v710  has  the  highest  sample  rate  of  the 
four  cameras  and  the  shortest  available  exposure  at  296  ns.  Recent  purchase  of  a  Phantom  v710 
for  the  DERF  ensures  its  availability  for  testing. 
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Table  A-1.  High-speed  camera  comparison 


Phantom  v7.0 

Phantom  vl2.1 

resolution 

Frame  rate 

Sample  rate 

Resolution 

Frame  rate 

Sample  rate 

Y 

Y 

X  (pix)  (pix) 

(frame/s) 

(pix/s) 

X  (pix)  (pix) 

(ffame/s) 

(pix/s) 

800 

600 

4796 

2.30E+09 

1280 

800 

6242 

6.39E+09 

640 

480 

7207 

2.21E+09 

1280 

720 

6933 

6.39E+09 

512 

512 

8213 

2.15E+09 

512 

512 

20978 

5.50E+09 

256 

256 

26143 

1.71E+09 

256 

256 

66997 

4.39E+09 

128 

128 

67796 

l.llE+09 

128 

128 

183250 

3.00E+09 

64 

64 

121212 

4.96E+08 

128 

64 

330469 

2.71E+09 

32 

32 

160000 

1.64E+08 

128 

8 

1000000 

1.02E+09 

Phantom  v710 

Fastcam  SA-5 

Resolution 

Frame  rate 

Sample  rate 

resolution 

frame  rate 

sample  rate 

Y 

Y 

X  (pix)  (pix) 

(ffame/s) 

(pix/s) 

X  (pix)  (pix) 

frame/s 

pix/s 

1280 

800 

7530 

7.71E+09 

1024 

1024 

1000 

1.05E+09 

512 

512 

25000 

6.55E+09 

832 

444 

20000 

7.39E+09 

256 

256 

79000 

5.18E+09 

512 

373 

50000 

9.55E+09 

128 

128 

215600 

3.53E+09 

256 

64 

300000 

4.92E+09 

128 

32 

685800 

2.81E+09 

128 

64 

420000 

3.44E+09 

128 

16 

1077500 

2.21E+09 

128 

24 

775000 

2.38E+09 

128 

8 

1400000 

1.43E+09 

64 

16 

1000000 

1.02E+09 

An  undesired  effect  of  detonation  is  the  light  generated  by  combustion.  Hydrogen/air 
emits  mostly  in  the  UV  band,  and  the  effect  on  imaging  is  small  (Fig.  A-7a).  Hydrocarbon  fuels 
emit  much  more  light  in  the  visible  range  (Fig.  A-7b),  and  the  light  obscures  the  cell  structure  of 
a  detonation. 
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Detonation 
front 

Self-luminance 

Weak  emission  from  H2-Air  Strong  emission  from  HC-Air  detonation 

Detonation 

Figure  A-7.  Light  emission  from  detonations 

Because  the  emission  depends  on  wavelength,  a  spectral  filter  can  block  the  light  from 
detonation  while  passing  light  from  another  source.  Figure  35  shows  the  emission  spectra  of 
hydrogen/air  and  acetylene/air  detonations.  Both  detonations  have  weak  emission  from  650  nm 
to  700  nm.  Meanwhile,  the  CMOS  sensors  employed  by  all  of  the  cameras  considered  are  near 
peak  sensitivity  at  700  nm  (Fig.  A-8c).  Filtering  for  such  a  narrow  band  of  wavelengths  will 
reduce  the  intensity  of  light  reaching  the  camera.  A  brighter  source  may  be  necessary  for 
sufficient  illumination. 
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H2+0. 502+1. 88N2  detonation  at  lOlkPa 


C2H2+Air  flame  at  lOlkPa 


Phantom  vl2.1  speetral  response:  blaek  line  denotes  monoehrome  model. 

Figure  A-8.  Spectral  emission  from  combustion  and  camera  sensitivity 

The  Phantom  v710  measures  position,  time,  and  speed  with  low  uncertainty.  The  bias 
uncertainty  in  position  is  a  function  of  the  image  resolution  and  exposure  time.  Again,  consider 
the  theoretical  stoichiometric  hydrogen/air  detonation  this  time  combined  with  the  Phantom  v710 
camera.  The  wave  speed  is  1971  m/s  and  the  cell  size  is  8.19  mm.  The  camera  settings  are  64 
by  64  pixel  resolution,  685800  frame/s,  and  296  ns  exposure.  Uncertainty  in  the  position  of  an 
object  is  the  root  sum  squared  of  the  uncertainties  due  to  pixel  size  and  the  distance  traveled 
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during  exposure.  The  uncertainty  due  to  pixel  size  varies  from  0.200  mm  to  0.0139  mm 
depending  on  the  selected  resolution  (32  pixels  and  128  pixels  respectively  over  the  200  mm 
wide  test  section).  The  total  uncertainty  ranges  from  0.0443  mm  to  0.204  mm.  The  pixel  size 
was  the  most  important  factor  when  setting  up  high  frame  rate  imaging  and  should  be  as  small  as 
allowable  for  sufficient  frame  rate.  The  bias  uncertainty  in  velocity  calculated  for  the  example  is 
15.5  m/s.  For  the  1100  m/s  to  2200  m/s  range  the  uncertainty  varies  from  7.48  m/s  to  18.1  m/s. 
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D-Series  Cases 

Case  D1 


Case  D3 
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R-series  Cases 
Case  R1 
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Case  R3 


Case  R4 


M-series  Cases 

Case  Ml 
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Case  M3 
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Case  M5 


Case  M6 
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Case  M7 
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Case  M9 


Case  MIO 
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Case  Mil 
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XI.  Appendix  C  -Video  Stills 


D-Series  Cases 

Case  D1 :  Straight  Channel  (0  =  0,  r  =  co) 

Run  1  Run  2  Run  3  Run  4 
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Case  D2:  0  =  90°,  r  =  2.0  mm 
Run  1  Run  2 


Run  3 


Run  4 
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Case  D2  (continued) 


Run  6 


Run  5 


Run  7 
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Case  D3:  0  =  15°  r  =  2.0  mm 

Run  1  Run  2  Run  3  Run  4 
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Case  D3  (continued) 
Run  5 


Run  6 


Run  7 


Run  8 
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Case  D3  (continued) 


Run  9 


Run  10 
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Case  D4:  0  =  90°,  r  =  25.4  mm 

Run  1  Run  2 


Run  3 


Run  4 
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Case  D4  (continued) 

Run  5 

Run  6 

Run  7 

Run  8 
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Case  D4  (continued) 


Run  9  Run  10 


Run  4 


XII.  R-Series  Cases 

Case  Rl:  p  =  45°,  xo  =  162.2  mm,  yo  =  -50.2  mm 

Run  1  Run  2  Run  3 
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Case  R1  (continued) 

Run  5 


Run  6 
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Case  R1  (continued) 


Run  9 


Run  10 
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Case  R2:  p  =  45°,  xo  =  84.7  mm,  yo  =  -50.8  mm 

Run  1  Run  2  Run  3 
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Run  4 
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Case  R2  (continued) 

Run  5 


Run  6 


Run? 


Runs 
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Case  R2  (continued) 


Run  9 


Run  10 
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Case  R3:  p  =  30°,  xo  =  80.1  mm,  yo  =  -50.8  mm 

Run  1  Run  2 


Run  3 


Run  4 
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Runs 


Case  R3  (continued) 

Run  5 


Run  6 


Run? 
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Case  R3  (continued) 


Run  9 


Run  10 
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Case  R4:  p  =  15°,  xo  =  0.0  mm,  yo  =  50.8  mm 

Run  1  Run  2  Run  3  Run  4 
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Case  R4  (continued) 


Run  5 


Run  6 
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Case  R5:  p  =  30°,  xo  =  80.1  mm,  yo  =  -50.8  mm,  r  =  25.4mm 

Run  1  Run  2  Run  3 


Run  4 
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Case  R5  (continued) 


Run  5 


Run  6 
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Case  R6:  p  =  -45°  xO  =  0.0  mm,  yO  =  50.8  mm 

Run  1  Run  2 


Run  3  Run  4 
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Case  R6  (continued) 

Run  5  Run  6  Run  7  Run  8 
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Case  R6  (continued) 


Run  9 


Run  10 
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Case  R7:  p=  -45°,  xO  =  43.0  mm,  yO  =  50.8  mm 

Run  1  Run  2  Run  3 


Run  4 
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Case  R7  (continued) 

Run  5 


Runs 


Run  6 


Run? 
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Case  R8:  p  =  -45' 

Run  1 


,  xo  =  169.3  mm,  yO  =  50.8  mm 

Run  2 


Run  3 


Run  4 
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Case  R8  (continued) 

Run  5 

Run  6 

Run  7 

Run  8 
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XIII.  M-Series  Cases 


Case  Ml 

Run  1 


Run  2 


Run  3 


Run  4 
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Case  M2 


Run  1 


Run  2 


Run  3 


Run  4 
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Case  M2  (continued) 

Run  5 


Run  6 


Run? 


Runs 
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Case  M2  (continued) 

Run  9 


Run  10 


Run  11 


Run  12 
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Case  M2  (continued) 

Run  13  Run  14 


Run  15 


Run  16 
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Case  M2  (continued) 


Run  18 


Run  19 
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Run  20 
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Case  M2  (continued) 


Run  2 1  Run  22 
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Case  M3 

Run  1  Run  2  Run  3  Run  4 
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Case  M3  (continued) 

Run  5 


Run  6 


Run? 


185 


Case  M4 

Run  1  Run  2  Run  3  Run  4 
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Case  M4  (continued) 

Run  5 

Run  6 

Run  7 

Run  8 
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Case  M5 

Run  1  Run  2  Run  3  Run  4 
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Case  M5  (continued) 

Run  5 


Run  6 


Run  7 


Run  8 
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Case  M5  (continued) 

Run  9 

Run  10 

Run  11 

Run  12 
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Case  M5  (continued) 
Run  13 


Run  14 


Run  15 


Run  16 
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Case  M6 

Run  1 


Run  2 


Run  3 


Run  4 
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Case  M6  (continued) 

Run  5 


Run  6 


Run  7 
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Case  M7 

Run  1 


Run  2 


Run  3 


Run  4 
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Case  M7  (continued) 

Run  5 


Run  6 


Run  7 


Run  8 
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Case  M7  (continued) 


Run  9 


Run  10 
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Case  M8 

Run  1  Run  2  Run  3  Run  4  Run  5 
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Case  M9 


Run  1 


Run  2 
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Case  M9  (continued) 

Run  3 


Run  4 
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Case  M9  (continued) 

Run  5 


Run  6 
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Case  M9  (continued) 

Run  7 


Run  8 
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Case  M9  (continued) 


Run  9 


202 


Case  MIO 


Run  1 


Run  2 
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Case  MIO  (continued) 

Run  3  Run  4 
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Case  MIO  (continued) 


Run  5 
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Case  Mil 


Run  1 


Run  2 
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Case  Mil  (continued) 

Run  3 


Run  4 
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Case  Mil  (continued) 

Run  5 


Run  6 
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Case  Mil  (continued) 

Run  7 


Run  8 
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XIV.  Appendix  D:  Source  Code 
Manual  pixel  selection  and  coordinate  output:  Markup2.m 

function  out  =  markup2(varargin) 

%  coords  =  markup2(franne) 

%  coords  =  markup2(Stack,  frameNo) 

%  coords  =  markup2(Stack,  frameNo,  options) 

% 

% 

%  markup2  allows  the  user  to  inteactively  select  pixels  from  an  image  and  returns 
%  the  coordinates  of  the  selected  pixels  as  a  Px2  matrix. 

%  Input:  frame(numeric)  -  a  matrix  of  pixel  values. 

%  Stack(lmStack)  -  an  image  stack 

%  frameNo(scalar)  -  frame  number 

%  options(ParamA/alue  pairs)  -  optional  arguments 
%  diff(logical)  -  In  diff  mode,  the  first  frame  of  the 
%  stack  is  subtracted  from  the  frame,  and  areas  outside 

%  the  region  if  interest  are  masked.  Diff  mode  does 

%  not  change  single  frame  inputs. 

%  cmap(string/nx3  numeric)  -  a  colormap  to  use  when  displaying 
%  intensity  images.  Does  not  change  true 

%  color  images. 

% 

%  Output:  coords(Px2  double)  -  a  list  of  coordinate  pairs  of  the  user  selected 
%  pixels 

%%  Input  checking  and  standardizing 

%number  of  arguments 
error(nargchk(1 ,6,nargin)); 

%  frame 
if  nargin==1 

if  ~isnumeric(varargin{1}) 
error('invalid  frame'); 
end 

Stack  =  lmStack(varargin{1}); 
frameNo  =  1 ; 
diffMode  =  false; 
cMap  =  'gray'; 
end 

%  Stack  and  frameNo 
if  nargin==2 

if  ~isa(varargin{1},'lmStack') 
error('invalid  ImStack'); 
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end 

Stack  =  varargin{1}; 

if  ~(isnumeric(varargin{2})  &&  isscalar(varargin{2})  && 
varargin{2}<=Stack.frameCount) 
error('invalid  frame  number'); 
end 

frameNo  =  varargin{2}; 
diffMode  =  false; 
cMap  =  'gray'; 
end 

%  frame  or  Stack  and  frameNo  with  parameters 
if  nargin>2 

if  isa(varargin{1},'lmStack')  &&  isnumeric(varargin{2})  &&  isscalar(varargin{2}) 
Stack  =  varargin{1}; 
frameNo  =  varargin{2}; 
diffMode  =  false; 
cMap  =  'gray'; 
argPtr  =  3; 

elseif  isnumeric(varargin{1}); 

Stack  =  lmStack(varargin{1}); 

frameNo  =  1 ; 

diffMode  =  false; 

cMap  =  'gray'; 

argPtr  =  2; 

else 

error('invalid  argument'); 
end 

%  loop  through  param/value  pairs 
for  iArg  =  argPtr:2:nargin 

%  switch  on  parameter  name 
switch  varargin{iArg} 

case  'difT  %  set  diff  mode 
try 

logical(varargin{iArg+1}); 
catch  ME 
clear  ME 

error('invalid  mode  argument'); 
end 
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diffMode  =  logical(varargin{iArg+1}); 

case  'cmap'  %set  colormap 
if  ~chkCMap(varargin{iArg+1}) 
error('invalid  colormap'); 
else 

cMap  =  varargin{iArg+1}; 
end 

%  reset  map  to  default  map:  jet 
if  strcmp(cMap, 'default') 
cMap  =  'jet'; 
end 

case  'clip'  %clip  color  scale 

assert(any(varargin{iArg+1}  ==  [0  1  false  true]),... 
'invalid  mode  argument'); 

otherwise 

error('invalid  parameter'); 
end 
end 
end 

%%  Apply  optional  differencing  and  color  map 

%  differencing 
if  diffMode  &&  frameNo>1 
frame  =  diff(Stack,frameNo); 
else 

frame  =  Stack(frameNo); 
end 

%  log  scale  and  clip  outliers 
%  frame  =  Iog10(frame-min(frame(:))+1); 

[n,b]  =  hist(frame(:),unique(frame)); 

cLim  =  [b(find(n==2,1, 'first')), b(find(n==2,1, 'last'))]; 

%  frame(frame<cLim(1))  =  cLim(1); 

%  frame(frame>cLim(2))  =  cLim(2); 

%  construct  colormap  and  convert  to  RGB 
if  strcmp(Stack.colorFmt, 'Monochrome'); 
cLen  =  length(unique(frame)); 
mn  =  min(frame(:)); 
mx  =  max(frame(:)); 

frame  =  round((cLen-1).*(frame-mn)./(mx-mn)+1); 
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map  =  eval([cMap,'(',num2str(cLen), ')']); 
frame  =  ind2rgb(frame,map); 


hmap  =  rgb2hsv(map); 

if  mean(hmap(:,2))<0.3  %low  saturation  (nearly  gray  scale) 
selColor  =  [1,0,0];  %red 

elseif  mean(hmap(:,3))  <  0.7  %low  value  (dark  colormap) 
selColor  =  [1,1,1];  %white 
else 

selColor  =  [0,0,0];  %black 
end 

else  %true  color  image 
selColor  =  [1,0,0];  %red 

end 

%%  Pre-Proc  frame  for  display  and  create  figure 
[n,m]  =  size(Stack); 

%  shade  outside  the  roi  (a  nice  soothing  blue) 

frame  =  cat(3,frame(:,:,1  ),frame(:,:,2),frame(:,:,3)+frame(:,:,3).*~Stack.roi); 
frame(frame>1)  =  1; 
cData  =  frame; 

%  create  custom  pointer 
cd  =  NaN(16); 
cd(8:9,:)  =  1; 
cd(:,8:9)  =  1; 
cd(8:9,6:11)  =  2; 
cd(6:11,8:9)  =  2; 
cd(7:10,7:10)  =  NaN; 


%  Initialize  figure 
hf  =  figure('lnterruptible','off'); 
ss  =  get(0,'ScreenSize'); 
op  =  [-7,33,ss(3)+16,ss(4)-24]; 
set(hf,'OuterPosition',op, . . . 

'Pointer', 'custom',... 

'PointerShapeCData',cd,... 

'PointerShapeHotSpot',[8,8]); 
hi  =  imshow(cData,map,'lnitialmagnification','fit'); 
hz  =  zoom(gcf); 
hp  =  pan(gcf); 
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axis  on 

title('Click  to  select  first  point  then  use  num  pad  to  select  more'); 

xlabel('Press  space  to  quit.'); 

pos  =  [.03, .03, .96, .96]; 

set(gca, 'Position', pos,'TickDir', 'in'); 

%%  Run  user  input  loop 
%  Solicit  first  point 
isDone  =  false; 
set(hz,'Enable','on'); 
waitfor(hz.  Enable, 'off'); 
while  -isDone 

k  =  waitforbuttonpressO; 
if  k 

out  =  []; 
return 

%  on  a  key  input  recycle 
else 

%  on  a  mouse  click  toggle  pan  and  zoom  modes  or  select  point 
if  strcmp(hz.Enable,'on') 
waitfor(hz,'Enable','off'); 

elseif  strcmp(hp. Enable, 'on') 
waitfor(hp,'Enable','off'); 

else 

curPt  =  get(gca,'CurrentPoint'); 
xLim  =  get(gca,'xLim'); 
yLim  =  get(gca,'yLim'); 
hWidth  =  round(diff(xLim)/2); 
hHeight  =  round(diff(yLim)/2); 
isDone  =  true; 
end 
end 

end  %while 

%  get  clicked  point 
j  =  round(min(max(curPt(1,1),1),m)); 
i  =  round(min(max(curPt(1,2),1),n)); 
points(1,:)  =  [i,j]; 
cData(i,j,:)  =  selColor; 

%  modify  figure 
set(hi,'cData',cData); 

%  add  labels 
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title('Use  numpad  to  add  pixels  to  selection.'); 
xlabel('Press  space  to  exit'); 

%  define  function  to  update  selection 
function  update() 
points(end+1,:)  =  [ij]; 
cData(ij,:)  =  selColor; 
set(hi,'cData',cData); 
end 

%  define  keypress  fen 
function  key_press(hf, event) 

%  runs  on  a  key  press  in  the  figure  should  update  cData  and  points  on  each 
%  key-press  and  exit  on  space 

%  interpret  key 

switch  event. Character 

case  '1'  %down-left 
j  =  round(nnax([1,j-1])); 
i  =  round(min(in,i+1])); 
updateO; 

case  '2'  %down 
%j=j 

i  =  round(nnin([n,i+1])); 
updateO; 

case  '3'  %down-right 
j  =  round(nnin([m,j+1])); 
i  =  round(min([n,i+1])); 
updateO; 

case  '4'  %left 
j  =  round(max([1,j-1])); 

%  i  =  i 
updateO; 

case  '5'  %  undo  last 
%reset  pixel 
cData(i,j,:)  =  frame(ij,:); 
set(hi,'cData',cData); 

%  remove  last  point  on  list 
points(end,:)  =  []; 
j  =  points(end,2); 
i  =  points(end,1); 
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case  '6'  %right 
j  =  round(min([m,j+1])); 

%  i  =  i 
updateO; 

case  7'  %up-left 
j  =  round(nnax([1,j-1])); 
i  =  round(max([1,i-1])); 
updateO; 

case  '8'  %up 
%j=j 

i  =  round(nnax([1,i-1])); 
updateO; 

case  '9'  %up-right 
j  =  round(nnin([m,j+1])); 
i  =  round(max([1,i-1])); 
updateO; 

case  ' '  %close  the  figure 
close(hf); 

end  %switch 

%  recenter  when  current  point  gets  close  to  edge 
i2  =  i-0.5; 
j2=j-0.5; 
nn2  =  m-0.5; 
n2  =  n-0.5; 

if  i2-yLim(1)<  10&&  i2>  10 
%  recenter  up 

yLim  =  [max([0.5  ,i2-hHeight]),... 

max([2*hHeight-0.5,i2+hHeight])]; 
set(gca, 'yLim', yLim) 
end 

if  yLim(2)-i2  <  10  &&  i2  <  n2-10 
%  recenter  down 

yLim  =  [min([n2-2*hHeight,i2-hHeight]),... 

min([i2+hHeight  ,n2  ])]; 

set(gca, 'yLim', yLim); 
end 

if  j2-xLim(1)  <  10  &&j2  >  10 
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%  recenter  left 

xLim  =  [max([0.5  ,j2-hWidth]),... 

max([2*hWidth-0.5J2+hWidth])]; 

set(gca,'xLim',xLim); 

end 

if  xLim(2)-j2  <  10  &&j2  <  m2-10 
%  recenter  right 

xLim  =  [min([m2-2*hWidthJ2-hWidth]),... 

min(02+hWidth  ,m2  ])]; 

set(gca, 'xLim', xLim); 
end 

end  %key_press 


set(hf,'KeyPressFcn',@  key_press); 
waitfor(hf); 

%  switch  from  ij  to  xy  ordering  (for  easy  plotting) 

out  =  [points(:,2),points(:,1)]; 

end 


% - 

%  Subfunctions 

% - 

function  tf  =  chkCMap(arg) 

%  validate  colormap  by  checking  against  the  list  of  built  in  maps  or  checking 
%  for  a  nx3  numeric  array 


if 

any(strcmp(arg, {'jet', 'hsv', 'hot', 'cool', 'spring', 'summer', 'autumn', 'winter', 'gray', 'bone', 'copp 
er', 'pink', 'lines', 'default','hilo'})); 
tf  =  true; 

elseif  isempty(arg) 
tf  =  true; 

elseif  isnumeric(arg)  &&  size(arg,2)==3 
tf  =  true; 
else 

tf  =  false; 
end 
end 
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Class  for  storing  and  manipulating  image  data:  ImStack 


classdef  ImStack 

%  ImStack  Create  a  multi-image  stack  object.  ImStack  tries  to  be  more 
%  useful  than  the  standard  arrays  when  working  with  images  and  movies. 

% 

%  Syntax: 

%  OBJ  =  ImStackO  creates  an  empty  image  stack 

%  OBJ  =  ImStack(array)  creates  a  stack  from  an  array 

%  OBJ  =  ImStack(fileName)  creates  a  stack  from  a  file 
%  OBJ  =  lmStack(fileName, frames)  creates  stack  and  loads  only  spac'd  frames 
% 


%  Input 

%  array  -  a  numeric  array  of  two,  three  or  four  dimensions: 

%  Height  -by-  Width 

%  Height  -by-  Width  -by-  Frames 

%  Heigth  -by-  Width  -by-  Colors  -by-  Frames 

%  fileName  -  a  string  containing  the  name  of  the  file  to  import.  Partial 
%  and  full  paths  are  also  accepted  as  long  as  the  file  is  on 

%  the  MATLAB  search  path. 

%  frames  -  a  numeric  vector  of  frame  numbers  to  load. 

% 

%  Output 

%  OBJ  -  an  image  stack  with  the  following  properties: 

%  height  -  image  heigth  in  pixels 

%  width  -  image  width  in  pixels 

%  frameCount  -  number  of  images  in  stack 
%  class  -  the  data  class  of  the  images  i.e.  double,  uintS,  etc. 

%  colorFmt  -  either  'Monochrome'  or  'RGB'  depending  on  the  format 
% 


% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 


and  methods: 

diff  -  returns  the  difference  between  each  image  and  the  first 
one  in  the  stack 

length  -  overloads  the  built-in  function  to  retrun  the  number  of 
frames 

norm  -  scales  the  images  from  0  to  1  converting  to  double  if 
necessary  (useful  for  the  imshow  function) 
read  -  load  images  from  file 

size  -  overloads  the  built-in  function  to  return  the  frame  size 
as  [width,  height,  frameCount]  for  monochrome  and 
[width,  height,  3,  frameCount]  for  RGB  images. 
convert_fmt  -  converts  the  imagesc  back  and  forth  between  monochrome 
and  RGB  formats.  R,  G,  and  B  channels  are  average  when 
converting  to  monochrome,  but  the  original  data  is  not  lost. 
Runnning  convert_fmt  again  restores  the  original  images. 
Converting  Monochrome  to  RGB  duplicates  the  original 


218 


%  imagesc  in  each  color  channel. 

%  This  method  is  useful  for  false  coloring,  and  uscaling. 

% 

%  A  note  on  indexing... 

% 

%  Retriving  properties  and  calling  methods  uses  the  standard  syntax,  but 
%  indexing  an  ImStack  object  returns  the  images  themselves  greatly  reducing  the 
%  complexity  of  code  needed  to  access  subsets  of  the  stack.  Indexing  works  as 
%  follows: 

%  Obj(scalar)  -  returns  frame  n 
%  Obj(vector)  -  returns  the  frames  in  the  vector 
%  Obj(array)  -  returns  a  subset  of  the  3D  [width, height,frame]  or  4D 
%  [width,  height,  color,  frame]  stack. 

% 

%  Chris  Stevens 
%  Last  Update:  11  Jul2012 

%  PUBLIC  PROPERTIES 
properties 
width 
height 
frameCount 
frameClass  =  'double' 
colorFmt  =  'Monochrome' 
roi 

medianPrame 

source 

bg 

times 

dt 

end 

%%  PRIVATE  PROPERTIES 

properties(Access  =  'private'.  Hidden) 
data 
dataPmt 
minVal 
maxVal 
end 

%%  PUBLIC  METHODS 
methods  (Access  =  'public') 


%%  ImStack  (constructor) 

function  This  =  ImStack(varargin) 
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%Use  existence  and  type  of  argument  to  determine  what  to  do 
switch  nargin 

case  0  %  empty  stack 
%  Initialize  properties 
This.frameClass  = 

This.colorFmt  = 

This. source  =  'Workspace'; 

case  1  %  array  or  file 
argin  =  varargin{1}; 
if  isnumeric(argln)  %  stack  from  array 
switch  ndims(argln) 

case  2  %single  mono  frame 
[This. height,... 

This.width]  =  size(argln); 

This.frameCount  =  1; 

This.frameClass  =  class(argln); 

This.dataFmt  =  'Monochrome'; 

case  3  %single  RGB  frame  or  multiple  mono  frames 
if  size(argln(3))  ==  3; 

%  can't  tell  from  the  argument  so  ask 

button  =  questdlg('3  Mono  frames  or  1  RGB  frame?',... 

'Color  Format:','Mono','RGB','Mono'); 

%  use  answer  to  set  properties 
switch  button 
case  'Mono' 

[This. height,... 

This.width,... 

This.frameCount]  =  size(argln); 

This.dataFmt  =  'Monochrome'; 
case  'RGB' 

[This. height,... 

This.width]  =  size(argln); 

This.frameCount  =  1 ; 

This.dataFmt  =  'RGB'; 
end 

else 

%  must  be  monochrome 
[This. height,... 

This.width,... 

This.frameCount]  =  size(argln); 

This.dataFmt  =  'Monochrome'; 
end 
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case  4  %true  color  frames 
[This. height,... 

This.width,... 

j  ■  ■  ■ 

This.frameCount]  =  size(argln); 
This.dataFmt  =  'RGB'; 


otherwise 

error('Could  not  creat  object,  invalid  array  dimension'); 
end 

This.frameClass  =  class(argln);  %  Use  same  class  as  input  array 
This.source  =  'Workspace';  %  All  arrays  are  sourced  from  the 

workspace 

This.data  =  argin;  %  Copy  input  array  to  "data"  property 

else 

%  stack  from  file(s) 

assert(ischar(argln)  ||  iscellstr(argln),'lnvalid  argument'); 

if  ischar(argln) 
fileList  =  {argin}; 
else 

fileList  =  argin; 
end 

%  check  for  existence 
for  iFile  =  1  :length(fileList); 

assert(exist(fileList{iFile},'file')==2, . . . 

'File  "%s"  not  found',... 
fileList{iFile}); 
end 

%  read  files 

This.source  =  fileList{1}; 

This  =  read(This); 

for  iFile  =  2:length(fileList) 

This  =  This.append(fileList{iFile}); 
end 

This.source  =  fileList; 
end 

%  File  name  with  a  frame  argument 
case  2 

assert(ischar(varargin{1})  &&  exist(varargin{1},'file')  ==  2,... 

'Bad  file  name  or  file  not  found'); 
if  isnumeric(varargin{2})  %  single  file  and  frame  argument 
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This.source  =  varargin{1}; 

This  =  This.read(varargin{2}); 
elseif  ischar(varargin{2})  %  two  files 
fileList  =  varargin; 

This.source  =  fileList{1}; 

This  =  read(This); 

for  iFile  =  2:length(fileList) 

This  =  This.append(fileList{iFile}); 
end 

This.source  =  fileList; 
else 

error('lnvlaid  argument'); 
end 

%  List  of  three  or  more  file  names 
otherwise 

for  i  =  1  margin 

assert(ischar(varargin{i})  &&  exist(varargin{i},'file') 
'Invalid  file  name'); 
end 

fileList  =  varargin; 

This.source  =  fileList{1}; 

This  =  This.readO; 

for  iFile  =  2:length(fileList) 

This  =  This.append(fileList{iFile}); 
end 

This.source  =  fileList; 
end  %switch 

%Set  min  and  max  properties 
This.minVal  =  min(This.data(:)); 

This.maxVal  =  max(This.data(:)); 

%  calc  the  region  of  interest  and  median 
This.roi  =  mask(This); 

This.medianFrame  =  median(This.data); 

This.bg  =  This.get_frames(1); 
end  %lmStack 


% - 

%  SIZE  -  Overload  built-in  SIZE  to  use  properties 

% - 

function  varargout  =  size(This) 
switch  This.colorFmt 
case  'Monochrome' 
if  nargout  ==  0 
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varargout{1}  =  [This. height, This.width,This.frameCount]; 
else 

varargout  =  {This. height, This.width,This.franneCount}; 
end 

case  'RGB' 
if  nargout  ==  0 

varargout{1}  =  [This.height,This.width,3,This.frameCount]; 
else 

varargout  =  {This.height,This.width,3,This.frameCount}; 
end 


end 

end 

o/o - 

%  LENGTH  -  Overload  built  ion  length  to  return  number  of  images  in  stack 

o/o - 

function  I  =  length(This) 

I  =  This.frameCount; 
end 


% - 

%  END  -  Overload  built-in  END  so  modified  indexing  works 

% - 

function  b  =  end(This,k,~) 
switch  k 
case  1 

b  =  This. length; 
case  2 

b  =  This.width; 
case  3 

switch  This.colorFmt 
case  'Monochrome' 
b  =  This. length; 
case  'RGB' 
b  =  3; 
end 
case  4 

b  =  This. length; 
end 

end  %end 


o/o - 

%  DIFF  -  returns  the  diffrence  between  the  spec'd  frame  and  the  first  frame 

o/o - 

function  frame  =  diff(This,frameNo, reference) 
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%  use  all  frames  if  not  set 
switch  nargin 
case  1 

frameNo  =  2:This.frameCount; 
reference  =  This.bg; 
case  2 

reference  =  This.bg; 


%  repeat  subtracted  frame  to  match  size 
fn  =  This.get_frames(frameNo); 
if  isscalar(frameNo) 
f1  =  reference; 
else 

switch  This.colorFmt 
case  'Monochrome' 

f1  =  repmat(reference,[1,1,numel(frameNo)]); 
case  'RGB' 

f1  =  repmat(reference,[1,1,1,numel(frameNo)]); 
end 
end 

%  subtraction  works  differently  on  uints  and  floats 
switch  This.frameClass 
%floating  point  subtraction 
case  {'double', 'single'} 
frame  =  fn-f1 ; 

%  uint  subtraction  scales  to  fit  within  range 
case  {'uint8','uint1 6','uint32','uint64'} 
ceil  =  intmax(This.frameClass); 
rawFrame  =  (fn/2+ceil/2-f1/2); 
floored  =  rawFrame-min(rawFrame(:)); 
scale  =  double(ceil)/double(max(floored(:))); 
frame  =  scale*floored; 
end 
end 


o/o - 

%  DIVIDE  -  Returns  the  specified  frame  divided  by  the  first  frame. 

o/o - 

function  frame  =  divide(This, frameNo) 

This.frameClass  =  'double'; 


%  return  all  frames  if  no  frameNo  given 
if  nargin  ==  1 
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frameNo  =  2:This.franneCount; 
end 

if  length(frameNo)  ==  1 

frame  =  This.get_frames(frameNo)./This.get_frames(1); 
else 

switch  This.colorFmt 
case  'Monochrome' 

repDims  =  [1,1,length(frameNo)]; 
case  'RGB' 

repDims  =  [1,1,1  ,length(frameNo)]; 
end 

frame  =  This.get_frames(frameNo)./repmat(This.get_frames(1 ), repDims); 
end 

frame(  isinf(frame)  |  isnan(frame))  =  0; 


O/o - 

%  NORM  -  returns  a  normalized  (0  to  1),  floating  point  version  of 
%  the  spec'd  frame. 

% - 

function  frame  =  norm(This,frameNo) 

if  ~strcmp(This.frameClass, {'double', 'single'}) 
This.frameClass  =  'double'; 
end 

%  get  frame(s)  to  norm 
This.frameClass  =  'double'; 
original  =  This.get_frames(frameNo); 

%scale  intensities 
mx  =  double(This.maxVal); 
mn  =  double(This.minVal); 
frame  =  (original-mn)/(mx-mn); 
end 


o/o - 

%  READ  -  import  frames  from  file 

o/o - 

function  This  =  read(This, frames) 

%  check  that  a  file  is  associated  with  the  object 
assert(~strcmp(This. source, 'Workspace'),... 

'%s  is  not  linked  to  a  file.  Cannot  read  frames',... 
inputname(l)); 
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%  Get  the  file  extension  and  the  list  of  compatible  file  types 
fName  =  This. source; 

[~,~,ext]  =  fileparts(fName); 

stillExt  =  imformats; 

videoExt  =  {'avi','mpg','wmv','asf  ,'asx'}; 

%  separate  actions  for  video  and  stills 
switch  ext(2:end); 
case  videoExt 
%  video  files 

File  =  VideoReader(fName); 
if  nargin  ==  1 
This.data  =  read(File); 
else 

This.data  =  read(File, frames); 
end 

case  [stillExt(:).ext] 

%  still  image  files 
This.data  =  imread(fName); 
otherwise 

error('Unsupported  file  type'); 
end 

%  set  height,  width,  frameCount,  and  colorFmt  properties 
if  ndims(This.data)  <=  3 
[This. height,... 

This.width,... 

This.frameCount]  =  size(This.data); 

This.dataFmt  =  'Monochrome'; 
else 

[This. height,... 

This.width,... 

j  ■  ■  ■ 

This.frameCount]  =  size(This.data); 

This.dataFmt  =  'RGB'; 
end 

%  update  frameClass 
This.frameClass  =  class(This.data); 
end  %read 


o/o - 

%  APPEND  -  append  frame(s)  from  other  sources 

o/o - 

function  This  =  append(This,fileName) 
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%  validate  file  name 

assert(exist(fileName,'file')==2,'File  not  found'); 

%  read  file 
oldData  =  This.data; 
oldClass  =  class(oldData); 
newData  =  imread(fileName); 
newClass  =  class(newData); 

%  convert  class  if  necessary 
if  ~isa(newData, oldClass) 
warning('lmStack:Append:rescaleOnAppend',... 

'bit  depth  mismatch,  interpolating  to  highest  bits/pixel'); 
switch  oldClass 

case  'double'  %  double/* 
newData  =  double(newData); 
case  'single' 
switch  newClass 

case  'double'  %  single/double 
oldData  =  double(oldData); 
otherwise  %  single/uint* 
newData  =  single(oldData); 
end 

otherwise 

switch  newClass 

case  'double'  %uint*/double 
oldData  =  double(oldData); 
case  'single'  %uint*/single 
oldData  =  single(oldData); 
otherwise  %uint*/uint* 

oldBits  =  str2double(oldClass(5:end)); 
newBits  =  str2double(newClass(5:end)); 
if  oldBits>newBits 

newData  =  cast(newData,oldClass).*(2''oldbits-1)/(2''newBits-1); 
elseif  oldBits<newBits 

oldData  =  cast(oldData,newClass).*(2''newBits-1)/(2''oldBits-1); 
end 

end  %  switch  newClass 
end  %  switch  oldClass 
end  %if 

%convert  color  format  if  needed 
switch  ndims(newData) 

case  2  %single  grayscale  frame' 
newFrameCount  =  1 ; 
newFmt  =  'Monochrome'; 
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case  3  %single  RGB  frame' 
if  size(newData,3)  ==  3; 

%  can't  tell  from  the  argument  so  ask 

button  =  questdlg('3  Mono  frames  or  1  RGB  frame?',... 

'Color  Format:', 'Mono', 'RGB', 'Mono'); 

%  use  answer  to  set  properties 
switch  button 
case  'Mono' 

newFrameCount  =  3; 
newFmt  =  'Monochrome'; 
case  'RGB' 

newFrameCount  =  1 ; 
newFmt  =  'RGB'; 
end 

else  %  must  be  monochrome 

newFrameCount  =  size(newData,3); 
newFmt  =  'Monochrome'; 
end 

case  4  %multiple  RGB  frames' 

newFrameCount  =  size(newData,4); 
newFmt  =  'RGB'; 
otherwise 

error('lnvalid  image(s)  in  file'); 
end  %switch 

%  default  to  RGB  if  formats  disagree 
if  ~strcmp(newFmt,This.dataFmt) 
warning('lmStack:Append:colorMismatch',... 

'New  color  format  does  not  match  old  format  defaulting  to  RGB'); 
if  strcmp(newFmt,'Monochrome') 
newData  =  repmat(newData,[1,1,3,1]); 
newFmt  =  'RGB'; 
else 

oldData  =  repmat(oldData,[1,1,3,1]); 
end 
end 

%concatenate  frames 
if  strcmp(newFmt,'RGB'); 

This.data  =  cat(4,oldData, newData); 
else 

This.data  =  cat(3,oldData, newData); 
end 
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%  update  propoerties 
This.frameClass  =  class(newData); 

This.dataFmt  =  newFmt; 

This.frameCount  =  This.frameCount+newFrameCount; 
end  %fcn  append 


o/o - 

%  MAX  -  overload  builtin  max  function 

o/o - 

function  mx  =  max(This) 

mx  =  cast(This.maxVal, This.frameClass); 
end 


% - 

%  MIN  -  overload  builtin  max  function 

% - 

function  mn  =  min(This) 

mn  =  cast(This.minVal, This.frameClass); 
end 


o/o - 

%  CONVERT_FMT  -  switch  between  mono  and  RGB  color  formats 

O/o - 

function  This  =  convert_fmt(This) 

%  change  the  type  to  the  opposite 
switch  This.colorFmt 
case  'Monochrome' 

This.colorFmt  =  'RGB'; 
case  'RGB' 

This.colorFmt  =  'Monochrome'; 
end 
end 


% - 

%  CONVERT_TYPE  -  change  the  output  class  for  functions  and  indexing 

% - 

function  This  =  convert_type(This,type) 

%  Change  the  class  of  indexed  output 

if  ~strcmp(type, {'double', 'single', 'uint8','uint1 6', 'uint32','uint64'}) 
error('Unsupported  data  class'); 
end 

This.frameClass  =  type; 
end 


o/o - 

%  SET_BG  -  sets  the  background  image  used  in  diff  and  div 
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function  This  =  set_bg(This,bgFrames) 

%  average  selected  frames  to  create  a  mean  background 
This.bg  =  mean(This.data(:,:,:,bgFrames),4); 

end 


end  %methods  (public) 

o/o - HIDDEN  PUBLIC  METHODS 

methods  (Access  =  'public',  Hidden) 


% - 

%  SUBSREF  -  Overload  normal  subscripting  to  return  frames  for  a  scalar  index 

% - 

function  b  =  subsref(This,s) 

%  SUBSREF  Implementing  the  following  syntax: 

%  obj() 

%obj(1) 

%obj([1,2,  3]) 

%  obj. property 
%  obj.method(args) 

switch  s(1).type 

case  '()'  %  Array  indexing 

b  =  This.get_frames(s(1).subs{1}); 
case '.' 

%  property  access 
switch  s(1).subs 
%  public 
case  'height' 
b  =  This. height; 
case  'width' 
b  =  This.width; 
case  'frameCount' 
b  =  This.frameCount; 
case  'frameClass' 
b  =  This.frameClass; 
case  'colorFmt' 
b  =  This. colorFmt; 
case  'source' 
b  =  This. source; 
case  'roi' 
b  =  This. roi; 
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case  'bg' 
b  =  This.bg; 

%  hidden 
case  'data' 
if  length(s)  >  1 

b  =  This.data(s(2).subs{:}); 
else 

b  =  This.data; 
end 

case  'min' 
b  =  This.min(); 
case  'max' 
b  =  This.max(); 
case  'medianFrame' 
b  =  This. medianFrame; 
case  'times' 
b  =  This.times; 
case  'dt' 
b  =  This.dt; 

%  method  access 
case  'difT 

if  strcmp(s(2).subs{1 
b  =  diff(This); 
else 

b  =  diff(This,s(2).subs{:}); 
end 

case  'divide' 

if  strcmp(s(2).subs{:},':') 

b  =  divide(This); 
else 

b  =  divide(This,s(2).subs{:}); 
end 

case  'norm' 

if  strcmp(s(2).subs{:},':') 

b  =  norm(This); 
else 

b  =  norm(This,s(2).subs{:}); 
end 

case  'read' 

b  =  read(This,s(2).subs{:}); 
case  'convert_fmt' 
b  =  convert_fmt(This); 
case  'convert_type' 

b  =  convert_type(This,s(2).subs{:}); 
case  'append' 
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b  =  append(This,s(2).subs{:}); 
case  'length' 
b  =  length(This); 
case  'size' 
b  =  size(This); 
case  'set_bg' 

b  =  set_bg(This,s(2).subs{:}); 
%  throw  controlled  errors 
otherwise 

if  numel(s)  ==  1 

error('Unknown  property'); 
else 

error('Unknown  method)'); 
end 

end  %switch  s(2) 
otherwise 

error('Syntax  error') 
end  %switch  s(1) 
end  %subsref 
end  %methods  (hidden) 

o/o - PRIVATE  METHODS - 

methods  (Access  =  'private') 


o/o - 

%  GET_FRAMES  -  Returns  frames  with  proper  class,  color  format,  and  subscripting 

o/o - 

function  outFrames  =  get_frames(This,args) 

%convert  color  format  if  needed 
fmt  =  strcmp('RGB',{This.dataFmt,This.colorFmt}); 
if  fmt(1  )==fmt(2)  %  same  format 

outFrames  =  This. data; 

elseif  fmt(2)  %  mono  data/rgb  frames 

temp  =  reshape(This.data,... 

[This.height,This.width,1,This.frameCount]); 
outFrames  =  repmat(temp,[1,1,3,1]); 
else  %  rgb  data/mono  frames 

outFrames  =  mean(This.data,3); 
end 

%convert  class 

outFrames  =  cast(outFrames,This.frameClass); 


%sub  sample  full  array 
rows  =  1  This. height; 
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cols  =  1  This.width; 
switch  This.colorFmt 
case  'RGB' 
colors  =1:3; 
isColor  =  true; 
case  'Monochrome' 
colors  =  []; 
isColor  =  false; 
end 

frames  =  1:This.frameCount; 
switch  length(args) 

case  1  %  single  full  frame 
frames  =  args(1); 
case  2  %  specified  pixels  only 
rows  =  args{1}; 
cols  =  args{2}; 
case  3  %  pixels  and  frames 
rows  =  args{1}; 
cols  =  args{2}; 
frames  =  args{3}; 
case  4  %fully  spec'd  true  color 
rows  =  args{1}; 
cols  =  args{2}; 
colors  =  args{3}; 
frames  =  args{4}; 
otherwise 

error('lndex  exceeds  dimensions') 
end 

if  isColor 

outFrames  =  outFrames(rows, cols, colors, frames); 
else 

outFrames  =  outFrames(rows,cols,frames); 
end 

end  %get_frames 
end  %methods 

end  %classdef 
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Function  to  mask  solid  objects  in  images:  roiMask 

function  roiMask  =  mask(This) 

%  roiMask  =  mask(This) 

% 

%  mask  returns  a  logical  array  which  is  true  within  the  region  of  interest  and 
%  false  elsewhere.  The  roi  is  the  light  region  of  the  first  frame  in  Stack  with 
%  the  largest  area. 

% 

%  Input:  Stack  -  an  image  stack  see  the  IMStack  class  for  more  info. 

%  Output:  roiMask  -  a  logical  array  the  same  size  as  the  frames  of  Stack  that 
%  is  true  within  the  roi  and  false  elsewhere. 

im  =  This.get_frames(1); 

%Use  intensity  thresh  to  separate  visible  areas  from  black 
raw_roi  =  im  >  max(im(:)).*0.1; 

%remove  undesired  sections 
roiMask  =  true(size(raw_roi)); 
s  =  regionprops(~raw_roi,'Area','PixelldxList'); 
area  =  cat(1,s.Area); 
pil  =  cat(1,s(area>200).  Pixel  IdxList); 
roiMask(pil)  =  false; 
end 

function  medFrame  =  median(data) 
m  =  (size(data,4)+1)*0.5; 
fs  =  sort(data,4);  %sort  frames 

if  mod(m,1 )  %  even  frame  count 

medFrame  =  (fs(:,:,:,m+0.5)+fs(:,:,:,m-0.5)).*0.5; 
else  %even  frame  count 
medFrame  =  fs(:,:,:,m); 
end 
end 
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Function  to  interpret  binary  video  data:  read_cine 

function  out  =  read_cine(varargin) 

%  out  =  read_cine(fileName) 

%  Reads  a  .cine  file  (Vision  Research  video  format) 

% 

%  Input: 

%  fileName  (string)  -  name  of  the  cine  file 

%  range  (2x1  numeric)  -  optional  range  of  frames  to  read  from  file 

%  option  (string)  -  Optionally  one  of  the  three  strings:  'ImStack',  'Array',  or  'Struct' 

%  which  specify  the  output  format.  ImStack  is  a  class  with 

%  some  rudimentary  analysis  methods.  Array  is  a  4D  array  of 

%  pixel  intensities.  Struct  is  a  structure  containing 

%  infromation  from  the  cine  file  as  well  as  pixel  values. 

%  If  not  specified,  read_cine  returns  a  stuct 

% 

%  Output: 

%  out  (varies)  -  The  struct  option  returns  with  the  following  fields: 

%  frameRate 

%  exposure 

%  frameCount 

%  version 

%  bitDepth 

%  width 

%  height 

%  colorFormat 

%  -  The  Array  option  is  double  class  and  4D  (height,  width,  color 

%  frame) 

%  Type  'help  ImStack'  for  information  about  the  class 

%%  Validate  arguments 
error(nargchk(1 ,3,nargin)); 
fileName  =  varargin{1}; 
switch  nargin 
case  3 

assert(any(numel(varargin{2})  ==  [0,2])  &&  isnumeric(varargin{2}), 'Invalid  range'); 
range  =  varargin{2}; 

assert(any(strcmpi(varargin{3},{'lmStack', 'Array', 'Struct'})), 'Invalid  option'); 
outclass  =  lower(varargin{3}); 
case  2  %  name  and  range  only 

assert(numel(varargin{2})  ==  2  &&  isnumeric(varargin{2}), 'Invalid  range'); 
outclass  =  'imstack'; 
case  1  %  name  only 
range  =  []; 

outclass  =  'imstack'; 
end 
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%%  Validate  and  open  file 

assert(exist(fileName,'file')  ==  2, 'File  not  found');  %  does  it  exist 
global  fid 

fid  =  fopen(fileName,'r'); 

assert(fid  ~=  -1  ,'File  unreadable');  %  did  it  open 

fileMarker  =  read(2,'CHAR'); 

assert(strcnnp(fileMarker,'Cr),'File  is  not  a  cine  file');  %  is  the  marker  correct 

%%  Read  cine  file  header 
headerSize  =  read('WORD'); 
compression  =  reacl('WORD'); 
version  =  read('WORD'); 
firstMovielmage  =  read('LONG'); 
totallmageCount  =  read('DWORD'); 
firstImageNo  =  read('LONG'); 
imageCount  =  read('DWORD'); 
offImageHeader  =  read('DWORD'); 
offSetup  =  read('DWORD'); 
offImageOffsets  =  read('DWORD'); 
triggerTime  =  read('TIME64'); 
fPos=  ftell(fid); 

assert(fPos  ==  headerSize, 'File  read  error:  headerSize  mismatch'); 

%%  Read  bit  map  info  header 

%  check  position  and  go  to  beginning  of  BITMAPINFOHEADER 
if  fPos  ~=  offImageHeader 
fseek(fid, OffImageHeader, 'bof); 
end 

biSize  =  read('DWORD'); 
biWidth  =  read('LONG'); 
biHeight  =  read('LONG'); 
biPlanes  =  read('WORD'); 
biBitCount  =  read('WORD'); 
biCompression  =  read('DWORD'); 
biSizelmage  =  read('DWORD'); 
biXPelsPerMeter  =  read('LONG'); 
biYPelsPerMeter  =  read('LONG'); 
biCIrUser  =  read('LONG'); 
biCIrImportant  =  read('DWORD'); 

%%  Read  setup  structure 
frameRate16  =  read('WORD'); 
shutter16  =  read('WORD'); 
postTrigger16  =  read('WORD'); 
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frameDelaylG  =  read('WORD'); 
aspectRatio  =  read('WORD'); 
contrast16  =  read('WORD');  %unused 
bright16  =  read('WORD');  %unused 
rotate16  =  read('BYTE');  %unused 
timeAnnotation  =  read('BYTE');  %unused 
trigCine  =  read('BYTE');  %unused 
trigFrame  =  read('BYTE'); 
shutterOn  =  read('BYTE');  %unused 

%  read  description  until  0x5343  ('ST') 
descriptionOld  =  read(2,'CHAR'); 
while  ~strcmp(descriptionOld(end-1  :end),'ST') 
descriptionOld  =  [descriptionOld, read('CHAR')]; 
end 

mark  =  descriptionOld(end-1:end); 
descriptionOld  =  descriptionOld(1  :end-2); 
length_  =  read('WORD'); 
binning  =  read('WORD'); 
sigOption  =  read('WORD'); 
binChannels  =  read('SHORT'); 
samplesPerlmage  =  read('BYTE'); 
binName  =  cell(8,1); 
for  i  =  1 :8 

binName{i}  =  read(1 1 , 'STRING'); 
read('BYTE'); 
end 

anaOption  =  read('WORD'); 
anaChannels  =  read('SHORT'); 
res6  =  read('BYTE'); 
anaBoard  =  read('BYTE'); 
chOption  =  read(8, 'SHORT'); 
anaGain  =  read(8, 'FLOAT'); 
anaUnit  =  cell(8,1); 
for  i  =  1 :8 

anaUnit{i}  =  read(5, 'STRING'); 
read('BYTE'); 
end 

anaName  =  cell(8,1); 
for  i  =  1 :8 

anaName{i}  =  read(10, 'STRING'); 
end 

iFirstImage  =  read('LONG'); 

dwImageCount  =  read('DWORD'); 

nOFactor  =  read('SHORT'); 

wCineFileType  =  read('WORD');  %#ok<*NASGU> 
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szCinePath  =  cell(4,1); 
for  i  =  1 :4 

szCinePath{i}  =  read(65, 'STRING'); 
end 

bMainsFreq  =  read('WORD');  %unused 

bTimeCode  =  read('BYTE');  %unused 

bPriority  =  read('BYTE');  %unused 

wLeapSecDY  =  read('WORD');  %unused 

dDelayTC  =  read('DOUBLE');  %unused 

dDelayPPS  =  read('DOUBLE');  %unused 

genBits  =  read('WORD');  %unused 

res1  =  read('INT');  %ignore 

res2  =  read('INT');  %ignore 

res3  =  read('INT');  %ignore 

imWidth  =  read('WORD'); 

imHeight  =  read('WORD'); 

edrShutter16  =  read('WORD'); 

serial  =  read('UINT'); 

saturation  =  read('INT'); 

res5  =  read('BYTE');  %ignore 

autoExposure  =  read('UINT'); 

bFlipH  =  read('BOOL'); 

bFlipV  =  read('BOOL'); 

grid  =  read('UINT'); 

frameRate  =  read('UINT'); 

shutter  =  read('UINT'); 

edrShutter  =  read('UINT'); 

post! rigger  =  read('UINT'); 

frameDelay  =  read('UINT'); 

bEnableColor  =  read('BOOL'); 

cameraVersion  =  read('UINT'); 

firmwareVersion  =  reacl('UINT'); 

softwareVersion  =  read('UINT'); 

recordingTimeZone  =  read('INT');  %reads  18000  should  be  -5 
cfa  =  read('UINT'); 

bright  =  read('INT')*10;  %  converted  to  sw  scale 

contrast  =  10'^(read('INT')/100);  %converted  to  sw  scale 

gamma  =  10''(read('INT')/100);  %  converted  to  sw  scale 

reserved  1  =  read('INT');  %ignore 

autoExpLevel  =  read('UINT'); 

autoExpSpeed  =  read('UINT'); 

autoExpRect  =  read('RECT'); 

wbGain{i}  =  read(4,'WBGAIN'); 

rotate  =  read('INT'); 

wbView  =  read('WBGAIN'); 

realBPP  =  read('UINT'); 
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convSMin  =  read('UINT'); 
convSMax  =  read('UINT'); 
filterCode  =  read  ('I  NT'); 
filterParam  =  read('INT'); 
uf  =  read('IMFILTER'); 
blackCalSVer  =  read('UINT'); 
whiteCalSVer  =  read('UINT'); 
grayCalSVer  =  read('UINT'); 
bStampTime  =  read('BOOL'); 
soundDest  =  read('UINT'); 
frpSteps  =  read('UINT'); 
frplmgNr=  read(16,'INT'); 
frpRate  =  read(16,'UINT'); 
frpExp  =  read(16,'UINT'); 
mcCnt  =  read('INT'); 
mcPercent  =  read(64, 'FLOAT'); 
ciCalib  =  read('UINT'); 
calibWidth  =  read('UINT'); 
calibHeight  =  read('UINT'); 
calibRate  =  read('UINT'); 
calibExp  =  read('UINT'); 
calibEDR  =  read('UINT'); 
calibTemp  =  read('UINT'); 
headSerial  =  read(4,'UINT'); 
rangeCode  =  read('UINT'); 
rangeSize  =  read('UINT'); 
decimation  =  read('UINT'); 
masterSerial  =  read('UINT'); 
sensor  =  read('UINT'); 
shutterNs  =  read('UINT'); 
edrShutterNs  =  read('UINT'); 
frameDelayNs  =  read('UINT'); 
imPosXAcq  =  read('UINT'); 
imPosYAcq  =  read('UINT'); 
imWidthAcq  =  read('UINT'); 
imHeightAcq  =  read('UINT'); 
description  =  read(4096, 'STRING'); 

%%  tagged  info  blocks 
fseek(fid,offSetup+length_,'bof'); 
if  (offSetup+length_)<offlmageOffsets 
readMore  =  true; 
while  readMore 

blockSize  =  read('DWORD'); 
type  =  read('WORD'); 
moreBlocks  =  read('WORD'); 
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switch  type 

case  1002  %Time  only  block 

imTimes  =  read(imageCount,TIME64'); 
case  1003  %Exposure  only  block 

ImExp  =  read(innageCount,'DWORD')./2'^32; 
case  1005  %  Binary  signal  block 

binSignal  =  read(blockSize{ctr}-8,'BTYE'); 
case  1006  %  Analog  signal  block 

anaSignal  =  read(blockSize{ctr}-8,'BYTE'); 
otherwise 

error('Undefined  block  type'); 
end 

readMore  =  moreBlocks; 
end 
end 

%%  pointers  to  images 
switch  version 
case  0 

pimage  =  read(imageCount, 'DWORD'); 
case  1 

pimage  =  read(imageCount,'INT64'); 
end 

%%  images 

%  Initialize  data  structure 
im  =  struct('annotationSize', {},... 

'annotation',!},... 

'imageSize',!},... 

'pixels',!}); 

%  Calculate  number  of  pixels  and  data  format 
if  biBitCount  >  16 

nPixels  =  3*biWidth*biHeight; 
classStr  =  'WORD'; 
elseif  biBitCount  >  8 

nPixels  =  biWidth*biHeight; 

ClassStr  =  'WORD'; 
else 

nPixels  =  biWidth*biHeight; 

ClassStr  =  'BYTE'; 
end 

%  Read  frames 
if  isempty(range) 

range  =  [1,imageCount]; 
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end 


if  range(1)<1 
range(1)  =  1; 
end 

if  range(2)>imageCount 
range(2)  =  imageCount; 
end 

nFrames  =  range(2)-range(1)+1 ; 
if  cfa  ==  0 

pixels  =  zeros(innHeight,imWidth,1 , nFrames); 
else 

pixels  =  zeros(imHeight,imWidth, 3, nFrames); 
end 

for  i  =  range(1):range(2) 
index  =  i-range(1)+1 ; 

im(index).annotationSize  =  read('DWORD'); 
forj  =  1:im(index).annotationSize-1*8 
im(index).annotation{j}  =  read('WORD'); 
end 

im(index).imageSize  =  read('DWORD');  %size  in  bytes  divide  by  two  for  DWORDS 
pixels(:,:,:, index)  =  flipud(reshape(read(nPixels,classStr),biWidth,biHeight)'); 
end 

%  Close  the  file 
fclose(fid); 

%  Build  output 
switch  outclass 
case  'imstack' 

out  =  ImStack(pixels); 
out.times  =  imTimes(:,6)-imTimes(1,6); 
out.dt  =  mean(diff(out.times)); 
out.source  =  fileName; 
case  'array' 
out  =  pixels; 
case  'struct' 

out  =  struct('fileName', fileName,... 

'frameCount', imageCount,... 

'frameRate',frameRate,... 

'exposure', shutterNs, . . . 

'edr',edrShutterNs,... 

'bitDepth',realBPP,... 

'height', imHeight,... 

'width', imWidth,... 
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'frames', pixels,... 

'cameraSerial', serial,... 
'triggerTime',triggerTime, . . . 
'imageTime',imTimes,... 
'imageExp',imExp); 
end 

end  %fcn 

%%  subfunctions 
function  out  =  read(varargin) 

%  read  the  file  and  return  data  of  a  certain  type 
error(nargchk(1 ,2,nargin)); 
global  fid 


if  nargin  ==  1 
count  =  1 ; 
type  =  varargin{1}; 
else 

count  =  varargin{1}; 
type  =  varargin{2}; 
end 

switch  type 
case  'BYTE' 

out  =  fread(fid,count,'ubit8'); 
case  'CHAR' 

out  =  fread(fid,count,'*char')'; 
case  'WORD' 

out  =  fread(fid,count,'ubit16'); 
case  'INT16' 

out  =  fread(fid,count,'*int16'); 
case  'SHORT' 

out  =  fread(fid,count,'int16'); 
case  'BOOL' 

out  =  logical(fread(fid,count,'ubit32')); 
case  'DWORD' 

out  =  fread(fid,count,'ubit32'); 
case  'UINT' 

out  =  fread(fid,count,'*uint32'); 
case  'LONG' 

out  =  fread(fid,count,'int32'); 
case  'I NT' 

out  =  fread(fid,count,'*int32'); 
case  'INT64' 

out  =  fread(fid,count,'*int64'); 
case  'FLOAT' 
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out  =  fread(fid, count, '*single'); 
case  'DOUBLE' 

out  =  fread(fid, count, '*double'); 
case  'STRING' 

out  =  fread(fid,count,'*char'); 
case  'TIME64' 

out  =  zeros(count,6); 
for  i  =  1  :count 

fraction  =  fread(fid,1  ,'ubit32')/2''32; 
seconds  =  fread(fid,1,'ubit32'); 

out(i,:)  =  datevec(seconds./(24*3600)+datenum('31  Dec  1969  18:00')); 
out(i,6)  =  out(i,6)+fraction; 
end 

case  'IMFILTER' 

out  =  struct('dim',[], 'shifts', [], 'bias', [],'coef',[]); 
for  i  =  1  :count 

out(i).dim  =  fread(fid,1,'int32'); 
out(i).shifts  =  fread(fid,1  ,'int32'); 
out(i).bias  =  fread(fid,1,'int32'); 
out(i).coef  =  fread(fid,25,'int32'); 
end 

case  'WBGAIN' 

out  =  struct('r',[],'b',[]); 
for  i  =  1  :count 

out(i).r  =  fread(fid,1 , 'single'); 
out(i).b  =  fread(fid,1 , 'single'); 
end 

case  'RECT' 

out  =  struct('r',[],'c',[],'h',[],'w',[]); 
for  i  =  1  :count 

out(i).r  =  fread(fid,1  ,'int32'); 
out(i).c  =  fread(fid,1  ,'int32'); 
out(i).h  =  fread(fid,1  ,'int32'); 
out(i).w  =  fread(fid,1  ,'int32'); 
end 
end 
end 
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Function  to  calculate  distance  between  two  curves:  cp_dist 

function  [dist,pos]  =  cp_dist(coords1,coords2,nnode) 

%  dist  =  cp_dist(curve1,curve2,mode) 

% 

%  cp_dist  runs  on  one  of  three  modes:  'closest',  'interp',  or  'combined' 

% 

%  In  'closest'  mode,  it  returns  a  vector  of  the  distance  from  curvel  to  curve2 
%  using  the  closest  point  on  curve2  to  each  point  on  curvel .  The  result  is 
%  the  same  length  as  curve  one. 

% 

%  In  'interp'  mode,  attempts  to  align  the  meaurement  better  to  the  normal  of  the 
%  curves.  Each  curve  is  interpolated  so  that  the  number  of  data  points  is  the 
%  average  length  of  the  two  curves,  then  the  distance  is  calculated  point  to  point 
% 

%  In  'combined'  mode  interpolates  each  curve  as  in  'interp'  mode  then  uses  the 
%  closest  interpreted  point  to  measure  distance. 

% 

%  In  both  cases,  cp_dist  returns  dist  and  pos.  Dist  is  a  nx3  array 
%  [x-distance,  y-distance,  magnitude].  Pos  is  a  2xn  array  of  the  origins  of  the 
%  distance  vectors. 

switch  mode 

case  {1, 'Closest', 'closest'} 

%  Closest  point 
n  =  length(coordsl); 
dist  =  zeros(n,3); 
for  iPt  =  1  :n 

d  =  sqrt((coords2(:,1  )-coords1  (iPt,1  )).'^2+(coords2(:,2)-coords1  (iPt,2)).'^2); 
[mag.ind]  =  min(d); 

dist(iPt,:)  =  [coords2(ind,1)-coords1(iPt,1),... 
coords2(ind  ,2)-coords1  (iPt,2), . . . 
mag]; 
end 

pos  =  coords  1 ; 

case  {2, 'Interp', 'interp'} 

%  Interpolated  point-to-point 
n1  =  length(coordsl); 
n2  =  length(coords2); 
nint  =  max([n1,n2]); 

c1i  =  [interp1q((1:n1)',coords1(:,1),linspace(1,n1,nlnt)'),... 
interpi  q((1  :n1  )',coords1  (:,2),linspace(1  ,n1  ,nlnt)')]; 


244 


c2i  =  [interp1q((1:n2)',coords2(:,1),linspace(1,n2,nlnt)'),... 
interpi  q((1  :n2)',coords2(:,2),linspace(1  ,n2,nlnt)')]; 

dist=  [c2i(:,1)-c1i(:,1),... 
c2i(:,2)-c1i(:,2),... 

sqrt((c2i(:,1  )-c1  i(:,1  ))/2+(c2i(:,2)-c1  i(:,2)).'^2)]; 


pos  =  c1i; 


case  {3, 'Combined', 'combined'} 

%  Interpolate  points  on  each  curve  then  use  closest  point  for  distance 
n1  =  length(coordsl); 
n2  =  length(coords2); 
nlnt  =  2*max([n1,n2]); 

c1i  =  [interp1q((1:n1)',coords1(:,1),linspace(1,n1,nlnt)'),... 
interpi  q((1  :n1  )',coords1  (:,2),linspace(1  ,n1  ,nlnt)')]; 

c2i  =  [interp1q((1:n2)',coords2(:,1),linspace(1,n2,nlnt)'),... 
interpi  q((1  :n2)',coords2(:,2),linspace(1  ,n2,nlnt)')]; 

dist  =  zeros(nlnt,3); 
for  iPt  =  1  mint 

d  =  sqrt((c2i(:,1)-c1i(iPt,1)).''2+(c2i(:,2)-c1i(iPt,2)).''2); 

[mag.ind]  =  min(d); 
dist(iPt,:)  =  [c2i(ind,1)-c1i(iPt,1),... 
c2i(ind,2)-c1i(iPt,2),... 
mag]; 
end 

pos  =  {c1i,c2i}; 
otherwise 

error('lnvalid  mode  argument'); 


end 
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Function  to  calculate  2D  polynomial  fitting  coefficients:  polyfitl 

function  C  =  polyfit2(x,y,z,  method) 

%  polyfit2  is  for  2-D  data  fitting  using  least  squares 
% 

%  USAGE:  C  =  polyfit2(X,Y,Z,  'method') 

%  where  an  output  vector  C  contains  the  bi-linear  or  bi-cubic 
%  coefficients  of  a  least-squares  polynomial  in  x  and  y,  and 
%  input  matrices  X,  Y,  Z  are  for  a  2D  function  z=f(x,y). 

% 

%  Here  'method'  can  be 
%  'linear'  -  bilinear  least  squares  fitting 

%  'cubic'  -  bicubic  least  squares  fitting 

%  n  -  binomial  of  order  n  least  squares  fitting 
%  Non-equally  spaced  (or  even  non-monotonic)  X  and  Y  are  permitted. 

% 

%  For  example,  generate  a  coarse  2D  curve  and  a  least  squares  fitting 

%  over  finer  mesh  (meshdom  with  Matlab  3.5  BUT  meshgrid  with  Matlab  4) 

%  X  =  0:10;  y  =  1 :9:  [x  y]  =  meshdom(x,y) ; 

%  z  =  sin(x.*y); 

%  xi  =  0:.25:10:  yi=2:.5:8  ;  [xi  yi]=meshdom(xi,yi): 

%  C  =  polyfit2(x,y,z,  'cubic'); 

%  zi  =  polyval2(C,  xi,yi,  'cubic'); 

%%  Check  arguments 
error(nargchk(2, 4, nargin, 'struct')); 

if  size(x)~=size(z), 

error('X  must  have  the  same  dimension  as  Z.'); 
end 

if  size(y)~=size(z), 

error('Y  must  have  the  same  dimension  as  Z.'); 
end 

if  ~ischar(method) 
n  =  method; 

method  =  [num2str(n),'th  order']; 
if  n  >=  numel(x) 

error('Order  must  be  less  than  number  of  data  points'); 
elseif  numel(n)  >  1 

error('Order  must  be  a  scalar'); 
end 
end 

%  Default  to  bilinear  fit 


246 


if  nargin<4 
method  =  'linear'; 
end 

%%  Calculate  coefficients 
X  =  x(:): 

y  =  y(:): 

z  =  z(:); 

len  =  length(z); 

%  Calculate  A  matrix 
switch  method 
case  'linear' 

A  =  [  ones(len,1 ),  x,  y,  x.*y  ] ; 
case  'cubic' 

A  =  [  ones(len,1),  x,  y,  x.*y,  x^2,  y.'^2,  (x.'^2).*y,  x.*(y.'^2),  xA3,  y.'^S] 
otherwise  %  nth  order  binomial 
n  =  n+1; 

A  =  zeros(length(x),sum(1  :n)); 


ctr  =  1 ; 
for  nx  =  0:n-1 
for  ny  =  0:n-1 
if  nx+ny  <  n 

A(:,ctr)  =  x.'^nx  .*  y/ny; 
ctr  =  ctr+1; 
end 
end 
end 
end 

C  =  A  \  z; 
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